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Abstract 


There are always differences between theoretical and experimental results in the study of lifting 
surfaces. Bounding box control volume measurements infrequently yield exact conservation of mass 
or consistent values for lift and drag coefficients. Numerically calculated wakes often differ from 
experimental data. Quite often, an empirical correction can be applied to fit theory to experiment 
to account for these differences. However, as the demands for state of the art foil design increase, 
fluid dynamicists are pressed to look carefully at these inconsistencies in order to improve current 
design and analysis methods. Using a Reynolds Averaged Navier Stokes (RANS) computer code and 
a highly refined fluid mesh, one can begin to explore the subtle characteristics of the fluid flow in 
the entire domain and the details of certain key regions around a foil. Specific areas of great interest 
are: flow around the trailing edge, flow within the boundary layer, wake profiles and the influence 
of tunnel wall boundaries in experimental facilities. 

The overall goal of this thesis is to resolve some of the discrepancies between theoretical results 
and experimental data. A computer code has been developed to generate the geometry for the fluid 
flow domain surrounding an arbitrary foil shape at a specified angle of attack in the MIT Marine 
Hydrodynamics Laboratory (MHL) water tunnel. This geometry is provided as input data for the 
RANS solver. A suite of software tools are developed to provide post processing analysis to compare 
the RANS solution with other numerical techniques and experimental measurements. 

Through the use of case studies, the numerical results of the RANS code are compared with 
recent MHL experimental data and other computational tools. A comparison is made between 
the experimental and RANS code results using a control volume analysis. Boundary layer and 
wake profiles are also compared. A correction scheme is developed to extrapolate experimental 
measurements to unbounded fluid flow. 
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Chapter 1 


Introduction 


1.1 Overview 


Recent application of advanced computational methods to design new hydrofoil! sec- 
tions for propeller blades has yielded inconsistent results. There is a need in the hy- 
drodynamic research community to undertake an in depth study of two-dimensional 
lifting surfaces since these sections are the fundamental building blocks of most three 
dimensional design methods. Additionally, a sound methodology should be applied 
to make effective comparisons between numerical and experimental results. Once the 
differences between numerical predictions, experimental measurements and real appli- 
cation are understood, those lessons can be incorporated into computational methods. 
Some improvement will come from better understanding of fluid flow within boundary 
layers near the foil surface. Other areas in need of refinement are the effects of wake 
diffusion and flow separation near the trailing edge. 

Ultimately, yielding these improvements will be a multi-step process that involves 
significant computational and experimental effort far beyond the scope of a single 
thesis. The work herein represents the first few steps towards the goal of developing 
a more complete understanding of the subtleties of fluid flow around two-dimensional 


foils. 


1 Throughout this thesis, the words hydrofoil, 2-D lifting surface and foil will be used as synonymous terms. 
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1.2 Background 


1.2.1 Why Is There Interest in 2-D Flow? 


Presently, the hydrodynamic research community has turned its attention back to 
the study of two-dimensional lifting surfaces. The renewed interest is derived from 
recent attempts to design advanced hydrofoil sections for propeller applications that 
demonstrate improved cavitation characteristics without incurring a significant drag 
penalty. 

Recent work in the MIT MHL[17] water tunnel by Jorde[16] and Kimball[19] on 
foils with non-traditional camber distributions and unique trailing edges indicate 
that the above goals may be attainable. However, they also point out that our 
knowledge of two-dimensional flow around foils is incomplete. Bloch[7] developed 
a methodology for adding corrections to inviscid propeller design codes to account 
for anti-singing trailing edges. Several issues have come to light. First, the current 
computer design codes work well for conventional foils. However, the current codes do 
not work well for foils with cupped or blunt trailing edges or for foils with advanced 


camber distributions designed to delay cavitation inception. 
1.2.2 A Recent Design Problem 


As part of an advanced technology demonstration(ATD), new high speed propellers 
were designed and tested at the Carderock Division, Naval Surface Warfare Center. 
The design effort required to achieve the improved propeller performance using ad- 
vanced blade sections was monumental. The design process involved iterating several 
times between using propeller design computer codes and model testing[3]. The costs 
and time associated with this type of design procedure are prohibitive. 

It is not practical to expect that any commercial operator buying a propeller 
for a ship could fund such a research and development effort. But, the commercial 


operators do want the advantages these new foil sections offer. The way to make 
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designing with advanced foil sections cost effective is to eliminate the costly iteration 
between design computer codes and model testing by developing a database that 
generalizes these new geometries as suggested by Bloch. Prudent application of model 
testing and RANS analysis for à broad geometry of foil shapes could provide the 
necessary data to formulate correction routines for standard propeller design codes 


without significantly increasing the computation time for a converged solution. 
1.3 Motivation for Modeling 2-D Foils With RANS 


An incompressible Reynolds Averaged Navier Stokes (RANS) solver was chosen for 
the bulk of the numerical analysis in this thesis. A RANS computer code solves 
the equations of motion for each fluid element throughout an entire domain. When 
compared to other solution methods, such as an inviscid panel method, RANS is very 
time consuming and requires large amounts of random access memory (RAM) and 
central processing unit (CPU) time. This is a significant disadvantage, especially 
when one wants to analyze several foils at several angles of attack and Reynolds 
numbers. But, there are many aspects of a RANS solution that are attractive. 

The voluminous output from a RANS computer code contains all of the flow char- 
acteristics for every fluid element throughout the entire domain such as: Velocity, 
Pressure, and Shear Stress. Provided the flow domain is discretized with sufficient 
resolution, you can also extract and measure subtleties of the flow like boundary layer 
velocity profiles and thickness, and flow separation under adverse pressure gradients. 
It is desirable to be able to study the characteristics of a viscous wake behind a hy- 
drofoil. Also, with the abundance of data available in the solution, it is easy to make 
a direct comparison between conditions calculated at a prescribed location in the flow 
field and those measured in a geometrically similar experiment. 

A large part of this study involves characterizing the differences between numerical 


solutions and experimental measurements. In an unbounded flow, such as a 2-D 
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foil fixed in a uniform stream, RANS codes provide a good validation for computer 
programs such as PAN2D [21] and XFOIL [10] which are both inviscid panel methods 
coupled with integral boundary layer solvers. Neither PAN2D nor XFOIL is currently 
capable of modeling a foil in a flow constrained by walls including the viscous effects 
of the boundary layers, which is the case in water tunnel experiments. A RANS 
code can serve as a liaison between unbounded numerical codes and the bounded 
case of experiments. For an equivalent foil geometry, the RANS code can be used to 
characterize the differences in lift, drag and other properties for both bounded and 


unbounded flows. A methodology for doing this is presented in Chapter 5. 
1.4 Objectives 


As was alluded to in Section 1.3, one overall goal in this thesis is to provide a scheme 
for feeding back results and measurements taken in experiments as corrections that 
can be applied to the computationally efficient inviscid panel methods. Several steps 
need to occur before this goal can be realized in a substantial way. Methods need 
to be developed to efficiently generate input files for the RANS solver. As will be 
shown in Chapter 3, accurate geometric representation of a flow domain is perhaps 
the most demanding part of obtaining the solution. Accordingly, a great deal of effort 
was devoted to ensuring that the RANS domain geometry was exactly the same as 
the inviscid computer codes and the experimental setups. Computer codes need to 
be developed to reduce the output data to a tractable and meaningful form. The 
following are the critical path or the enabling objectives used to achieve the above 


goals: 


1. Review recent MIT MHL Water Tunnel experimental results for advanced section 


two-dimensional hydrofoils. 


2. Write a computer code that generates the fluid flow geometry input files for the 


RANS solver. 
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3. Analyze foils in unbounded flows using the RANS code and PAN2D. 
4. Analyze foils in bounded flows using the equivalent geometry of the MIT MHL. 


9. Develop corrective factors to apply to the inviscid solvers which are based on the 


differences found between the RANS and experimental results. 
6. Demonstrate how RANS solutions can be used to formulate experimental test plans. 
1.4.1 Rapid Generation of Flow Geometry 


To support ongoing research in the MIT Marine Hydrodynamics Laboratory, a two- 
dimensional lifting surface analysis tool is developed in this thesis to accurately model 
the fluid flow around hydrofoils in the water tunnel. It is very useful to have good 
predictions of a hydrofoil’s performance characteristics prior to conducting an exper- 
iment. Thorough empirical evaluation of hydrofoils requires taking many measure- 
ments with varying geometry and Reynolds numbers. It is common to take measure- 
ments for a single foil geometry at several angles of attack and Reynolds number. 
Therefore, one requirement for any computational tool used in the MHL is that it 
be easy to change the foil angle of attack and other test conditions such as scale and 
Reynolds number. Part of the work in this thesis is the development of the computer 
code FIT2D(Foil In Tunnel, Two-Dimensional). FIT2D is an interactive fluid mesh 
geometry generator. The user can arbitrarily specify: angle of attack, grid resolu- 
tion, and Reynolds number as well as many other lesser parameters. The output files 
generated are used as input for a Poisson equation (V? = Const) grid refinement 


program. After the grid is refined, ıt becomes input for a RANS solver for analysis. 
1.4.2 Resolving Differences Between Experiments and Numerics 


In sections 1.3 and 1.2.1, it is asserted that to make the current foil design computer 
codes work better, corrections should be incorporated to account for the physical 


effects that are not modeled in computationally efficient inviscid solutions. Kimball 
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found that there is uneven boundary layer growth on tunnel walls due to the presence 
of a foil that is generating lift[19]. This difference between the upper and lower walls 
affects the flow and the resulting lift and drag measurements. So there are really 
two sides to the correction scheme. One is using RANS to quantify the effects that 
the tunnel walls have on the foil lift and drag measurements as a result of uneven 
boundary layer growth and the potential flow imaging effect. The other is using the 
experimental measurements to apply corrections to computer codes for physics that 
are not captured by the numerics such as re-attachment of separated flows, transition 


back to a laminar boundary layer, and vortex shedding phenomena. 
1.5 Organization 


In Chapter Two, all of the relevant theory of two-dimensional lifting surfaces is out- 
lined. It provides the mathematical statement of the lifting problem. Similarly, 
Chapter Three is a detailed description of the numerical aspects of the computa- 
tional solution. 

Chapter Four outlines the methodology for post processing all of the RANS data. 
It provides an overview of the suite of programs that were developed and used to 
calculate lift and drag, and analyze wake profiles and boundary layers. 

Chapter Five compares and contrasts the RANS output to other numerical meth- 
ods and experimental measurements. Results of case studies for two different foils are 
presented. Figures, Graphs, and Tables are used to demonstrate the differences and 
similarities that occur. A methodology is proposed for applying corrections hydrofoil 
computer programs to gain better agreement with precise experimental results and 
real applications. A demonstration of how RANS solutions can be used to formulate 
experimental test plans is presented. 

Chapter Six is a summary of results and conclusions. Future topics and directions 


for research in this area are discussed. 
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Chapter 2 


Theory 


2.1 Overview 


This study focuses on 2-D fluid flow, therefore it is appropriate that this chapter 
outlines and develops the pertinent fluid mechanics theories for lifting surfaces and 
boundary layers. These theories have been previously derived in many texts and 
technical papers. The developments and derivations contained in this thesis represent 
the specific application to the problem of thin non-cavitating hydrofoils in bounded 
and unbounded incompressible two-dimensional flow domains. 

First, the basic lifting problem is defined as a boundary value problem and the 
geometry is specified. Subsequent sections draw attention to aspects of lifting surface 
theory that are fundamental to the computational methods employed and to the 


physical phenomena that are modeled. 
2.2 Two-Dimensional Lifting Flows 


Lifting surfaces have many applications in marine hydrodynamics. Two dimensional 
foil sections are used extensively in the design of rudders, keels, skegs, propellers, and 
other appendages, such as ducts, where a desired force is generated normal to the 
onset flow. 

This study is limited to a “typical” lifting surface as described in Newman [24] 


where the thickness is much smaller than the chord length (t(z)«C ). Further, only 
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incompressible two-dimensional flows are considered where the span of the foil is of 
infinite length and the resultant flow over the surface is chord-wise in direction. Figure 
2-1 shows how the foil geometry 1s defined with respect to the coordinate axes and 


the onset fluid flow. The foil is fixed in a stationary reference frame. The velocity 





Figure 2-1: 2-D Foil Coordinate System 


field (Uinf) propagates in the direction of the positive z-axis, and accordingly the 
leading edge of the foil points opposite to the onset flow. Although the pivot point is 
arbitrary for an unbounded flow, it becomes important when walls are present close 
to the foil surface. In this work, care was taken to specify the exact x,y location of 
the pivot point for comparing numerical results to those obtained in the MHL water 
tunnel. The angle of attack (a) is defined as the angle that an imaginary line drawn 
from the leading edge to the trailing edge makes with the z-axis. Here it is assumed 
that the “nose to tail” line intersects the fore and aft end of the mean camber line 
at the zero angle of attack. If the trailing edge is blunt, beveled or has some other 
treatment such as a splitter plate, the mean camber line ends at the midpoint of a 
line drawn from the upper and lower points of the trailing edge surface. The mean 


camber line is defined as: 


um(2) = 5(ya(2) + y(z)) (2.1) 
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Foil thickness is defined as: 
yılz) = yulz) — yılz) (2.2) 
Steady state flow is assumed for all cases. The Kutta Condition[18] applies at the 


trailing edge requiring that the fluid velocity be finite. With the above limitations 


and the stated geometry, the following boundary value problem results: 


NE 0 (2.3) 

V — 0,on y,(z) and у(х) (2.4) 

Уф « co,at the trailing edge, Kutta Condition (2.5) 
V = 0,on walls if present (2.6) 

Уф = Ving i, as T,y >00 (2.7) 


Equations 2.3 - 2.7 are the complete mathematical statement of the lifting problem for 
a stationary foil with onset flow. Two types of solutions to this problem are pursued. 
A Reynolds Averaged Navier Stokes solver is used to solve the steady state viscous 
flow around the foil. The other method is a potential flow vortex lattice method 
where the no-slip rigid boundary conditions of Equations 2.4 and 2.6 are relaxed. In 


inviscid potential flow, the wall and foil surface boundary conditions become: 
V-n=0 (2.8) 


Now that the boundary value problem is defined for the viscous and inviscid problems, 


the solutions can be formulated. 
2.3. Linearized Two-Dimensional Theory 


The use of linear potential theory is very powerful for analyzing lifting flows. A simple 
example of a lifting potential flow is a point vortex fixed in a uniform stream as shown 


in Figure 2-2. The linear potential for a free stream with a point vortex located at 
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SU. ds: 

ф- Пит + Im>-(log(z 0.) (2.9) 
Ihis potential is a solution to the Laplace equation (2.3) throughout the flow field 
except at the location of the point vortex. Using the principle of superposition, a 
more complicated flow can be formulated by adding a group of simple flows together 


that satisfy the kinematic boundary conditions and Laplace. 





Figure 2-2: Point Vortex in Uniform Stream 


2.3.1 Lift 


In potential flow the lifting force acting on a foil is obtained using the Kutta-Joukowsky 
theorem, which states that for any two-dimensional body, moving with constant ve- 
locity in an unbounded inviscid fluid, the hydrodynamic pressure force is directed 
normal to the velocity vector and is equal to the product of the fluid density, body 
velocity and the circulation about the body[24]. The mathematical statement of this 


theorem applied to thin foils is: 
L = $ (p-p..)da = pU Pu dar pUT (2.10) 


In a potential flow that is formulated by a distribution of discrete point vortices, 


application of 2.10 ıs a simple matter of summing up the strengths of the point 
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vortices enclosed by a simple contour. 
2.3.2 Drag 


In a pure inviscid potential flow solution, there is no pressure or form drag. However, 
when fluid viscosity is added to the solution, a viscous boundary layer(section 2.7) is 
formed near the foil surface. Drag on the foil is present as a result of two phenomena. 
First, shear stresses in the boundary layer cause skin friction drag. Second, due to the 
presence of the boundary layer, the pressure recovery at the trailing edge is incomplete 
which causes a drag force. The pressure drag is hard to calculate or measure because 
the pressure changes that cause the drag are isolated to a small area near the trailing 
edge of the foil. Also when compared to the lift, drag is usually a very small quantity. 
The errors associated with computing the drag are nearly equal to the drag itself. 
Two other methods of computing the drag are presented later in sections 2.5.1 and 
2.5.2. These methods are not as accurate as direct integration of pressure and shear 


stress but provide additional verification of the integrated quantities. 
2.3.3 Lift and Drag by Pressure Integration 


When the viscous effects are included, Kutta-Joukowsky does not apply for calculating 
lift. The drag is no longer zero. However, if the pressure and shear stress in the fluid 
adjacent to the surface(S) are known, the lift and drag can be obtained by direct 
pressure integration on the foil surface around the perimeter of the foil. The lift and 


drag forces on the foil due to pressure integration are: 
p / Р-п,45 (9.11) 
S 


and 


Ер = | Р.пуйѕ (2.12) 


where S is the foil surface and P is the fluid pressure at the foil surface. In addition 


to the normal pressure to the surface, there is a viscous shearing stress (Tas) acting 
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tangent to the surface. In a similar fashion to the above equations, shearing forces can 
be calculated by integrating around the foil and resolving into z and y components. 
Typically, the y component of the shearing forces is neglected because the integral 
quantity of these forces is insignificant compared to the y component of the pressure 


forces. The drag force on a foil due to viscous shearing forces is: 
boa / Tre’ ted (2.13) 
The total drag on the foil is the sum of the pressure and viscous drag terms, or: 
Ер = Ер, + Ер, (2.14) 
2.3.4 Viscous Effects on Lift 


A detailed overview of the viscous effects on lift is contained in Coney[9]. Some 
pertinent aspects of his discussion are presented here to provide needed clarity. As 
shown in Figure 2-3’, a foil in the presence of a real fluid flow will have a thin viscous 


boundary layer formed on its surface. This, in essence, alters the effective thickness 


Displ t 
Suction Side са 


Pressure Side 


Figure 2-3: Viscous Boundary Layer on 2-D Foil 


and camber distribution of the foil. The camber is reduced on the aft portion of 
the foil. In general, for the same foil geometry, the difference between the potential 
flow and viscous flow solution is that the change in camber distribution will cause an 
apparent reduction in lift coefficient. Boundary layer characteristics are covered in 


more detail in section 2.7 
! Figure 2-3 created by Scott D. Black 
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2.4 Vortex Lattice Methods(VLM) 


The basic concepts of vortex lattice theory are contained in Newman [24] and Ker- 
win [18]. In light of the long computation time to obtain a RANS solution, it is a 
prudent step to incorporate a vortex lattice analysis code into the computer code that 
generates the geometry for the RANS solver. This serves two purposes. First, it pro- 
vides an initial estimate for the lift coefficient(C;) and the circulation distribution(y) 
on a foil. Second, it provides the dividing streamline in the wake by extraction of 
the field point velocity due to the circulation and free stream potential. This is de- 
sirable because the numerics of the RANS solver works better if the fluid domain 
discretization is adapted to the location and trajectory of the wake.(See Section 3.7) 

Vortex lattice methods are simple and efficient. The earliest vortex lattice method 
is attributed to Falkner[12] as referencedby Kerwin[18]. The derivation presented here 
is analogous to Kerwin's method. The primary difference is that the actual mean line 
of the foil is used vice the linearized foil surface. 

The formulation of the vortex lattice method 1s basic. The mean camber line of 
a given foil is determined using equation 2.1. The camber line is broken up into 
a discrete number of elements of panels. The panels are spaced using a “cosine” 
method where the panels at the leading and trailing edge of the foil are smaller than 
the middle panels. A point vortex of unknown strength I’, is located at the mid-point 
of each of the m panels. Control points where the kinematic boundary conditions are 
imposed are located at the end of each panel. Figure 2-4 is an example of a vortex 


lattice distributed on a mean camber line. 
2.4.1 Cosine Spacing 


Cosine spacing is commonly used in many lifting surface and vortex lattice methods 


as the spacing algorithm for distributing point vortices and control points. It is a 
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Figure 2-4: Vortex Lattice on a 2-D Mean Line 


simple method where an auxiliary variable 3 1s defined such that: 


De 


2 (1 — cos $) (2.15) 





о) «== 


where 3 — Q is the leading edge of the mean line and 3 — 7 is the trailing edge. In 
the auxiliary coordinate system, the interval of the arc is divided into N equal panels 
each of $ — r/N length. To establish the locations for the point vortices and the 


control points the following equations are used: 








sin) S == [1 — с =, (2.16) 
san) — [1 — ссз(— (ЭЛТ) 


2.4.2 Obtaining the Vortex Lattice Solution 


The boundary condition on the foil is prescribed by equation 2.8. The strength of 
the point vortices are to be determined such that this boundary condition is satisfied 


at each of the control points. First it is instructive to have the expression for the 
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velocity field induced at a point z,y by a point vortex located at the origin[24]. 


Г и-]т 
This equation can easily be converted to the velocity field induced by a point vortex 
located at z,, y, on a control point located at £e, Ye by making a simple substitution 


of the difference in the coordinates for x and y. 


ІП (Ye Еу уу) — j(t. E no) 
МЕ c due to v — — 2n [(z. e z,)? i: (Ye ДА "H (2-19) 


lo extend the expression to apply to the vortices and control points distributed as 
per equations 2.16 and 2.17, simply insert the appropriate indices. Then the solution 
of the vortex lattice can be formulated into a mathematical statement. The velocity 
induced at the n™ control point by all the N vortices dotted with the normal vector 


is equal to mznus the free stream velocity dotted with the normal vector. 


Ip (п) tm) = H(ee(n)=2(m)) до роден 
a "ld am ula m) mn (2.20) 


Equation 2.20 can be written for each of the control points. This represents a set 
of N simultaneous equations which can be solved for the unknown vortex strengths 
(Ta). Separating out the [,, terms as a separate component, the terms that are left 
formulate a square matrix that represents the induced velocity caused by each point 
vortex on every control point. This matrix is called the influence coefficient matrix 


(A). Equation 2.20 can be cast in matrix form as: 
A-T=U (D 


where T is the array of point vortex strengths and U is the array of boundary condi- 
tions at the control points. The A matrix can be reduced using Gaussian elimination 
or any convenient algorithm for matrix inversion[4]. Once the matrix is inverted, the 


vortex strengths can be determined. 
DE ERU (2.22) 
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This completes the treatment of the 2-D vortex lattice method for a foil in an 
unbounded fluid. The next aspect to consider is the development of wall corrections 


for when a vortex lattice is located near a wall. 
2.4.3 Method of Images Applied to VLM 


Many of the cases analyzed involve 2-D foils which are in close proximity to walls. 
This is the case for testing a foil in a water tunnel. Therefore, the unbounded vortex 
lattice method needs to be modified to account for the presence of the walls. There 
are two options for modeling the walls. The first is to put vortex lattice panels on 
the walls in a similar fashion to panelizing the foil in the unbounded method. The 
second is to use the method of images where the walls are treated as “mirrors” and 
imaginary image foils are placed outside the flow дотаіп [28]. 

Adding panels to the walls 1s not an attractive option. Adequate representation of 
the walls with vortex panels requires many more panels than are already used for the 
foil. This adds significantly to the computer code solution time. Furthermore, the 
boundary condition at the walls (V -n = 0) is only satisfied at the individual control 
points, rather than continuously along the boundary. Another issue with paneling 
the walls is deciding how far up and down stream to extend the panelization scheme. 
The required distribution of wall vortex panels varies and is dependent upon many 
parameters. There is no obvious generalized scheme for implementing a rule based 
algorithm to determine the minimum required number of panels and their respective 
spacing on the walls. 

A better way to represent the walls in an inviscid flow solver is to use the method 
of images when the geometry permits. In this case, the foil is located between two 
parallel walls. Therefore the image geometry is very simple. Using the wall as a 
symmetry plane, an image body is placed on the opposite side of the wall from the 
actual body. For a vortex situated near one wall the image formulation is simple - - 


add one image vortex[28]. However with two parallel walls, the image geometry is a 
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little more complicated. The two parallel walls create an infinite number of reflective 
planes. This geometry, with the first pair of vortex lattice images, is shown in Figure 


2-5. Newman[24] provides a treatment for the similar case of a potential flow point 
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Figure 2-5: Vortex Lattice Near Walls With Images 


source situated between two walls located at y = +5b. Implementing this for the case 
of a distribution of point vortices is analogous. This requires the addition of an infinite 
array of image vortices at y = +6, +26, +30,---, to satisfy the boundary condition at 
the walls. A closed form solution can be formulated for the infinite array, but this 
Is unnecessary since each subsequent image pair has rapidly diminishing impact on 
the solution. For all practical purposes, a converged solution can be obtained with a 
small number of image pairs. This is demonstrated in section 3.5.2. The important 
thing to remember with this solution scheme is that the image vortices are the exact 
same strength as the original vortices. Therefore, no additional unknowns are added 
to equation 2.21. The only additional burden in computation is the addition of the 


influence of the vortex images to the influence coefficient matrix (A). Lastly, using 
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the image representation eliminates the problem of deciding how far up and down 
stream to extend the walls. Using images, the wall boundary condition is exactly 


satisfied for —oo « z « оо. 
2.5 Momentum Theory for Lift and Drag Calculations 


It has been shown that it 1s a simple task to calculate the the lift and drag by direct 
integration on the foil surface when using numerical foil analysis tools. However, in 
water tunnel experiments this is very difficult to do. It is very time consuming and the 
level of complexity of instrumentation required to take the surface measurements is 
impractical. The normal methods for measuring lift and drag for foils in the MIT MHL 
are “bounding box” velocity measurements using laser Doppler velocimetry(LDV)[17]. 
This is an application of fluid momentum theory. Bounding box calculations are used 
in this thesis to make comparisons between numerical calculations and experimental 
measurements. Bounding box contours can also be extracted from RANS solutions 


as a check of the surface pressure integration results. 
2.5.1 Bounding Box Analysis 


Bounding box methods are a way of calculating the forces acting on 2-D hydrofoils 
expressed in terms of integrals of the velocity flow field along a closed contour sur- 
rounding the foil[20]. Figure 2-6 is an example of a closed integration contour around 
a foil. Integrating around the contour C with an outward pointing unit normal vector 


n, the momentum flux M, out of C 1s: 
MS pp а(а:п)9з (2:23) 


where q is the fluid flux across C and p 1s the fluid density. The momentum flux 
through the surface of the foil is zero, since it 1s by definition a rigid boundary. The 
momentum flux of the volume of the fluid V enclosed on the outside by C and on 


the inside by the foil surface is given by equation 2.23. Using Newton's 3rd law, this 


33 








Figure 2-6: Rectangular Contour Around 2-D Hydrofoil 


momentum must be balanced by an opposing force F which is given by: 
piss $ Pnds + $ rtds — F; — Fp (2.24) 


When taking measurements around foils in the MIT MHL with the LDV apparatus, 
a rectangular contour is frequently chosen. This geometry is also convenient for 
extracting data from a RANS solver or from vortex lattice field point calculations. 
The contour is chosen such that it is sufficiently far from the foil so the flow is 
considered inviscid?. Therefore, all the contributions due to the shear terms (7) are 


zero and the pressure terms (P) are evaluated using the Bernoulli equation|20]. 
2.5.2 Drag by Wake Defect 


Although in principle, the direct application of equations 2.23 and 2.24 should yield 
adequate results, this is not always the case for drag measurements. À special treat- 
ment must be applied in order to get a better estimate for drag by bounding box 


methods[20]. First, some ancillary variables should be defined for the z-component 


? A minor exception to this is the small region where the wake passes through the downstream side of the box. 
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of the velocity in different regions of the flow along the right hand side of the contour: 


VEN, above the wake 


u“ = u + Au, in the viscous region of the wake 


-2 
М 
сл 
— 


DNE below the wake 2 


The first order calculation for drag is expressed as: 


Fp = pU ins I. 4 EE — u)dy (2.26) 


This can be further simplified using equation 2.25, resulting in: 


EDU (Au)dy (227) 


accross wake 
This equation is correct up to the first order of Au, as long as the contour which 
crosses the wake is sufficiently down stream such that Au < Ujn; holds. In the MIT 
MHL, the optical limits of the tunnel view port window do not allow for this to happen 
when large foils(c > 30cm) are tested. Therefore, a second order correction must be 
applied to equation 2.27. Typical measurements for bounding boxes around large 
foils have wake defect velocities which range as low as 0.3U;,; to 0.5U;,;. Kinnas[20] 
has shown, using equation 2.24 as a starting point, that the drag including the second 


order correction is: 


Шр = рб | (Au)dy — p (Auy dy (2.28) 


accross wake accross wake 
These equations for bounding boxes apply to both unbounded and bounded flows. 
Direct surface integration and bounding box analyses almost never agree exactly. In 
experiments there is the inherent uncertainty associated with taking a measurement. 
Additionally, the flow is usually complicated by some unsteady phenomena that oc- 
curs such as vortex shedding. In numerical computations there are usually small 
differences associated with the how the calculation scheme is implemented. When 


comparing experiments and computer results, the lift measurements usually agree 
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quite closely. Drag results tend to vary widely. This may well be due to the fact that 
drag, in comparison to all other flow characteristics, is a small quantity. Accordingly, 


it is a hard quantity to measure with great precision and consistency. 


2.6 Reynolds Averaged Navier Stokes Equations 


2.6.1 Derivation of the RANS Equations 


One cannot discuss the RANS equations without first reviewing the equations from 
which they originate. The Navier-Stokes Equations are the equations of motion for 
a viscous fluid. The derivation can be found in many texts, of which, Sabersky[28] 
provides a clear step by step formulation. The N-S equations stem from Newton’s 
law of motion and Newton’s law of viscous friction. The equations apply to viscous 
compressible flow with varying viscosity. The N-S equations also apply to turbulent 


flow. In vector form the equations are: 


= = - УР +vV’v+ т/7Ө = (2.29) 
The Mach number of the cases studied is sufficiently low such that the fluid is incom- 
pressible. Therefore the © term is zero. 

Considering the case of 2-D turbulent flow, the velocity of the fluid can be expressed 


as a time averaged quantity with a small fluctuating term added to it. Similarly, the 


pressure can be expressed as a time averaged pressure with a small fluctuating term. 
и›и=й+и 
о= 0+0 
PERET (2.30) 


The fluctuating terms u’,v’ and P’ are defined such that their time averaged value 
is zero. Substituting equations 2.30 into the N-S equations, yields(written for x 


component, similar for y): 


Ou ди да? т „дии " ди? E ди? " дау” * Ou'd " Qu'v' 
Ot s Ot st Ox Ox Ox Oy Oy Oy Oy 
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= 108-108 q yen (С "дейи, 
pOr  pór П ае “атэ dy?” 





(2.31) 


If the time averaged is taken for this equation, and the continuity equation is applied 


{ог о апа v the equation simplifies to: 


Du 19Р 19 да -- m ди Ny! 232 
Dt рдт pdr Cogo S p Oy onu du 


This equation is exactly equation 2.29 with the addition of two terms pu”? and pu'v”. 
The time averaged pressure and velocity satisfy the original N-S equation. In laminar 
flow these terms are zero. These extra terms are important in turbulent flow. They 


are commonly called the turbulent stresses or Reynolds stresses. 
2.6.2 Turbulence Modeling 


With the addition of the Reynolds Stress terms to the N-S equations, more equations 
must be included in the solution to evaluate the extra unknown terms. As these terms 
represent the turbulent characteristics of the flow, the extra equations required are 
usually called turbulence models. The word model is used because the formulation 
of the equations is typically based on empirical studies of turbulent flow. There are 
several models available for use in solving the RANS equations. Three of the most 


common derived models are: 


e Prandtl mixing length “zero equation” model 
e The “two equation” k — € model 


e Baldwin-Lomax algebraic model 


All the above methods employ the concept of eddy viscosity vr. The models use 
empirical relations to establish a value for vr which is fed into the RANS equations. 
The models used by the incompressible RANS solver in this thesis are the k — e and 
the Baldwin-Lomax Model. 

The Prandtl model uses a single parameter to determine the eddy viscosity. It 


is called the mixing length /. The mixing length is the characteristic length of the 
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mean eddy size which is much smaller than the fluids mean free path. It can also 
be thought of as the length normal to the main flow where momentum exchange is 


occurring. The Prandtl mixing model takes form: 


Dots (=) (2.33) 


This equation directly relates the eddy viscosity to the mean flow. Therefore, no 
additional unknown terms are added to the problem formulation. The mixing length 
is usually determined by the experimental results for flows which are similar physically 
to those being modeled numerically. 

A derivation of the k — e model is presented in White[30]. This model is based on 
dissipation. This requires the addition of two equations to the equations governing 
the fluid flow. One is a turbulent kinetic energy equation and the other is the tur- 
bulent dissipation equation. Included in these equations are five empirically derived 
constants. These constants can be varied so that they can be applied to different 
types of flows such as wakes and jets. If the solution does not predict the flow rea- 
sonably, then the model and its parameters probably do not match the physics of the 
chosen flow|29]. The model is not intended for flows with high curvature. In lifting 
surface analysis, there are situations where the flow has a high curvature near the 
leading edge or un-separated flow adhering to an anti-singing trailing edge. For this 
reason, the model was not extensively used in this thesis. 

The Baldwin-Lomax algebraic model is very attractive because of its computa- 
tional efficiency and accuracy[29, 30]. The model is appropriate for bodies at modest 
angles of attack where flow separation is minimal or is not expected to occur. This 


eddy viscosity model defines an effective viscosity ve. 
Ve = V + VT (2.34) 


The flow is divided into two regions. Different equations for the calculation of the 


turbulent eddy viscosity vr apply in each region. In the inner region the Prandtl 
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turbulent model(equation 2.33) is used with the van Driest[30] damping model for 
computing the mixing length. In the outer region a different set of equations applies. 


The eddy viscosity is calculated by[29]: 
ppm у Еке) (2.35) 


In the equation, K, Cop, Fake, and Fxıcs(y) are a set of constants and coefficients 
calculated from the local mean flow characteristics and the distance from a boundary. 
Equations 2.32(written for both u and v), 2.33, and 2.35 formulate a set of four 
equations with four unknown quantities. They are the velocities u and v, the pressure 


P and the eddy viscosity vr. 


2.7 Boundary Layer Flows 


2.7.1 Characterization of Boundary Layers 


A boundary layer, as the name implies, is a thin layer of fluid adjacent to a boundary 
such as a foil surface or a wall. Within the boundary layer the fluid undergoes a 
transition from the free stream velocity outside the boundary layer to at rest on the 
boundary. The majority of the viscous effects are confined to boundary layers. 

To employ a RANS solver effectively, it is important to discretize a flow domain 
so that an appropriate number of fluid elements are contained within the boundary 
layer to capture the essential viscous effects. In modeling flow around lifting surfaces, 
failure to adequately capture the flow in the boundary layer typically results in an 
overestimation of lift coefficient and underestimation of drag. If the boundary layer 
is missed completely the results tend toward the inviscid potential flow results. Since 
capturing the boundary layer flow is so important to the accuracy of the RANS 
solution, it is important to be able to estimate the boundary layer characteristics 
beforehand. The initial estimates for the boundary layer dimensions are based on 
flat plate theory which are also suitable for application to lifting surfaces operating 


at moderate angles of attack. 
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Many texts provide relevant discussion of boundary layer theory and character- 
istics, of which Newman[24] and White[30] are well suited to the level of discussion 


required here. 
2.7.2 Laminar Boundary Layer Flows 


Consider the case of a flat plate of length / in a two-dimensional fluid domain with 


uniform flow U. Equation 2.29 reduces to the following: 


u eue 

Ox Oy 
Ou Ou 10P ди Hu 
“Oe "Dy раз "дз + ду 


паши СИЕ Aff de 
dx Oy | рдт От Oye 


(2.36) 
The boundary conditions that apply for this flow are that u — U just outside the 
boundary layer and then u — 0 at the surface through a small distance ó. Within 


this small distance ó above the plate, the following characteristics are apparent: 
д 
5. =0(>) 
Oy д 


ди U 
эт = б\т) 


which implies: = > = 230) 


With u and v being 0 on the plate and in light of the relationships in 2.37, 2.36 
reduces to the following equations: 


Qu Qu 10P ди 
ИО = 


д2 "dy ps ду 
LEO (2.38) 

р ду 
The following conclusions can be drawn from equation 2.38: First, the pressure change 
across the distance 6 is not significant. Second, the pressure gradient in the x-direction 


inside the boundary layer is the same as the fluid outside[24]. In the absence of a 


pressure gradient in the z-direction, these equations can be solved numerically to 
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yield the Blasius solution for a flat plate. The concept of boundary layer thickness 
is introduced as the distance ó where the velocity u — 0.99U. In this formulation, 
the boundary layer thickness grows with the one-half power of the local Reynolds 


number(fe,). The local coordinate z is referenced from the leading edge of the plate. 


1 
ie (т) (2.39) 
U 
Another quantity that can be computed is the displacement thickness or 6*. This 
defines an effective thickening of the body due to viscous effects which corresponds 


to a lost quantity of fluid flux within the boundary layer. 
Бе) | (1 = z) dy г 1.72 (=) ў 2.40 
(ту p, dv 7 (2.40) 


Associated with the д^ flux reduction, there is a corresponding loss in fluid momentum 


and a momentum thickness which is calculated by: 


"m | 5 (1 » =) dy & 0.664 G (2.41) 


To calculate the drag on the flat plate due to viscous forces, the shear stress(T;,) at 


the plate must be determined. 


Qy 


Once the shear stress 1s known along the plate, the total frictional drag can be com- 


Tes (5) = 0.332pU? Rex? (2.42) 
y=0 


puted using equation 2.13. 

All of the above results for a flat plate can be extended to a general 2-D body 
without loss of validity, provided that the radius of curvature is much larger than the 
boundary layer thickness. 

The pressure gradient in the z-direction does affect the boundary layer. The 
boundary layer will become separated from the body at a point when the following 
conditions are met: the shear stress is zero; upstream of that point the tangential 
velocity is positive; and downstream of that point it is negative. A separated region 
is characterized by the streamlines breaking away from the body downstream. A 


separated wake is an area of relatively high vorticity and low pressure. 
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2.7.3 Turbulent Boundary Layers 


If the Reynolds number is sufficiently high, the laminar boundary layer will transition 
or change to a turbulent boundary layer at some distance aft of the leading edge of the 
foil. In the laminar boundary layer the velocity profile is smooth and regular and pre- 
dictable. At some point, the flow becomes disturbed and unstable. Transition occurs. 
After this, the boundary layer will typically remain turbulent along the remainder of 
the surface. The boundary layer remains relatively thin. The exact transition point 
can only be estimated and its location is dependent on Reynolds number and local 
roughness of the surface. Typically in experimental practice, turbulent transition is 
forced at a prescribed location near the leading edge of the foil. This is normally 
accomplished by the addition of surface irregularities such as raised bumps at the 
desired point of transition. Near the surface, inside the boundary layer, there is a 
viscous sub-layer that is small compared to the overall boundary layer thickness. The 
region between the viscous sub-layer and the outer fluid is called the turbulent core. 

The fluid flow in a turbulent boundary layer is more complicated than the laminar 
case. Mixing occurs at the interface between the free stream and the boundary 
layer. The basic momentum equations that apply to the laminar case apply to the 
turbulent boundary layer as well. At the surface the shear stress boundary condition 
still applies(t,, x Velocity gradient). The flow in the turbulent core and at the 
interface with the outer fluid is complicated. The assumption that the slope of the 
velocity profile provides an adequate value for shear stress(equation 2.42) is no longer 
valid[28]. As a result, empirical derivations of turbulent shear stress near a wall 
are typically used. One of the most common formulations is the 1/7'* power law 


approximation: 


1/7 
ые (=) (2.43) 


Ur и 


After making the substitutions u, — ,/T/p, u = Uing at y = 4, the following equations 


42 





result for shear stress at the surface(7,,) and turbulent boundary layer thickness(ó): 





TW 
ша EN : | 
т, = 0.0463 ( =) 002, (2.44) 
1/5 4/5 
§ = 0.058 | — (=) | 
| =) - (2.45) 


where the term a is determined by the computation of the integral: 


| и и y 
= а = 2 


These equations are included as they were used as a check for validating the quantities 





determined from the RANS results. À more complete discussion of turbulent layers 


and the development of the associated equations is contained in Sabersky[28]. 


2.7.4 The y* Parameter 


Finite element cell spacing near any no-slip boundary needs to be much smaller than 
the boundary layer thickness 6. This is necessary to accurately capture the velocity 
profile completely through the boundary layer to the no-slip surface. The distance 
that the close-in grid points are measured from the no-slip surface are measured in 


multiples of y* where: 


M 
| 
SIS 


(2.47) 


Ur = я (2.48) 


Based on Anderson[2], Black[6] and present RANS methodology in the MIT-MHL, 


and 


y+ values for the first three cells off the no-slip surface should be « 1, € 2, and < 4 
respectively. Additionally, there should be another 10 to 15 cells spaced from y* z 4 
to y* zz 100. It should be noted that failure to adhere to this tight spacing criteria 


can result in large errors. 
2.7.5 The Law of the Wall 


The Law of the Wall is a commonly used formulation for determining the velocity pro- 


file characteristics near a wall[28]. The main assertion is that the velocity profile(u) 
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is determined by the conditions at the маЩ таз), Ше fluid properties(p, v) and the 
distance from the wall(y). Using the Buckingham Pi theorem, two dimensionless 
parameters can be formulated from these quantities|23]. One is a non-dimensional 
velocity and the other a non-dimensional distance. There is an implied simple func- 
tional relationship between the two parameters: 


= ЕК (=) (2.49) 


и 


As a matter of convenience, the terms in equation 2.49 are often replaced by u* = 
f(y*). Several approximate formulas have been deduced based on experimental data. 
Of those available, the Spalding formula is quite suitable because it fits experimental 
data well from y* zz 100 all the way down to the surface[30]. 


(ки?) (кит)? 


E : (2.50) 


= T 


In the above equation, values for the constants « and B are 0.41 and 5.0 respectively. 


This formula is plotted in Figure 2-7. The Spalding formula for the Law of the Wall is 


1 
log,, (Y^) 


Figure 2-7: Spalding Formula for Law of the Wall 


used for validation of RANS calculations in selected cases to ensure that the velocity 
profiles close to the wall are correctly captured. One measure of the quality of a 


numerical finite element grid is how well the calculated boundary layer profiles agree 
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with the Spalding formula. It can be concluded that if the grid spacing conforms to 
the y* spacing criteria and the computed boundary layer profiles conform to the Law 


of the Wall, then the numerical grid is an adequate model for the flow. 
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Chapter 3 


Numerical Methods Development 


3.1 Overview 


Several 2-D foil analysis computer codes are employed in this study. The majority of 
the analysis was conducted using the RANS flow solver DTNS2D[29]. Other codes 
used are PAN2D[22] and XFOIL[11]. 

DTNS2D is a generalized 2-D domain incompressible RANS fluid flow solver devel- 
oped at the U.S. Naval Surface Warfare Center, Carderock Division. DTNS2D uses a 
finite volume finite difference formulation. In a finite difference solution method, the 
flow domain is divided into a set of discrete control volumes. Within each element, the 
governing equations of the flow domain are solved in their integral form. It is capable 
of modeling foils in bounded and unbounded domains. It should be noted that there 
are several RANS computer codes that could be applied to the cases studied here. 
However, DT NS2D was a readily available code within the MIT-MHL that required 
no modification for employment in this thesis. It has been used previously for ana- 
lyzing lifting surfaces with results that agree with experimental data[25]. Therefore, 
it was deemed adequate for this study as well. 

PAN2D and XFOIL use potential flow panel methods for calculating the inviscid 
flow characteristics around a foil in an unbounded 2-D domain. The solution is then 
coupled with an integral boundary layer method to determine the viscous properties 


of the flow[15]. The two panel methods were used to validate DTNS2D results for 
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unbounded flows. The two panel methods are currently not capable of modeling foils 
bounded by tunnel walls. 

In this chapter the development process for the FIT2D computer code that provides 
the fluid geometry input files for the DTNS2D solver is presented. This involves 
discretizing the domain into a mesh of four sided polygons finite volume elements. 
Grid generation schemes and methodology are discussed. As a subsection of FIT2D, a 
2-D vortex lattice computer code[18] was implemented to obtain an inviscid solution 
for the lift coefficient of a foil in unbounded and bounded flows. The vortex lattice 
solution is also used to grow a dividing streamline downstream from the trailing edge 


to define a RANS zonal grid boundary. 


3.2 Overview of Grid Generation 


3.2.1 Fundamental Concepts 


The overlying governing directive of grid generation is simple: discretize the domain 
into a sufficient number of elements such that the essential physics of the flow is 
captured. The implementation of this 1s difficult. Figure 3-1 is a representation of a 
discretized fluid domain surrounding a hydrofoil inside a tunnel. It is important to 
notice that the grid spacing is not uniform. The size of the individual elements must 
be varied depending on the expected local characteristics of the flow. For example, 
near a wall or on the surface of the foil, the vertical spacing must be quite small 
in order to accurately capture the velocity gradient within the boundary layer. For 
hydrofoils it is also important to cluster elements near the leading edge to capture the 
strong pressure gradients associated with the flow stagnation point. At the trailing 
edge, clustering is used where separated flow is expected. Additionally, the wake 
behind the foil has steep velocity gradients that slowly dissipate with the flow. The 
clustering scheme is advantageous if the small elements are located within the wake 


region. Conversely, there are regions in the flow where very little is occurring. Velocity 
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Figure 3-1: Sample of the discretized flow domain of a foil in a water tunnel. (Only 1/3 of the total 
grid lines are reproduced in the lower frame) 


and pressure gradients are near zero. In these areas, there is no reason to have finely 
spaced cells. 

In Figure 3-1, the element dimensions vary throughout the domain. The change 
in size from largest to smallest element is approximately four orders of magnitude. 
This sample representation requires 41,238 elements. In a grid representation of an 
unbounded flow domain around a foil, the element dimensions can range as much a 
six orders of magnitude. Of course one could simply find the smallest dimensional 
requirement for an element (eg. foil leading edge) and set all elements to that size and 


capture the essential physics of the flow. The problem is, for the same sample above, 


48 








the grid would require 6x10!? elements! In a typical numerical solver, the solution 
time Г 15 proportional to the number of elements squared (T' x n?). For the same 
results, the processor time to obtain a solution has grown by a factor of 10!?. Clearly 
this 1s not a practical solution. So, there are two conflicting requirements which must 
be met in order to produce results which are both practical and accurate. The first is 
to maximize the number of elements to capture all the pertinent physics. The other 
is minimize the number of elements to obtain the minimum solution time yet also 
obtain correct results. To satisfy both of these driving forces, the element sizes must 
vary throughout the domain in order to capture the local effects and minimize the 


overall number of elements. 
3.2.2 Geometric Limitations of DTNS2D 


Before proceeding with the description of the FIT2D program, it is important to spec- 
ify the geometric limitations for the DTNS2D input files. In the DTNS2D program 
the individual elements are grouped into zones. There can be any number of zones in 
DTNS2D. Figure 3-2 is a typical zone for DTNS2D. À zone consists of four sides. The 
sides can consist of any shape to enclose the zone so long as the four sides together 
comprise a simple connected region. À simple closed region has the property that an 
arbitrary closed curve lying in the region can be shrunk continuously to a point in 
the region without passing outside of the region|14]. Each pair of opposite sides must 
have the same number of elements along the side. Adjacent sides are not required to 
have equal number of elements. The elements within each zone cannot overlap any 
other elements. Every element within the zone must have four sides of finite non-zero 
length. Each zone side has a defined boundary condition (eg. no-slip, tangent, free 
stream, continuity of velocity gradient or pressure gradient). The boundary condition 
along a side can not change. Adjacent sides can have different boundary conditions. 
An example of a DTNS2D zone model is presented in Figure 3-3. This type of model 


is frequently called an “H” grid. 
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0.7 Typical Wall Boundary 
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Figure 3-2: Sample geometry of a DTNS2D zone. (Only 1/3 of the total grid lines are reproduced) 


3.3 Development of the Computer Code FIT2D 


3.3.1 Functional Requirements 


Recent research projects in the MIT-MHL water tunnel have demonstrated a re- 
quirement to use a RANS computer code to model flow around hydrofoils[19]. Some 
experiments on foils with unique camber distributions yielded interesting results. It 
is anticipated that the RANS output could help provide some insight into the results 
that were being observed. 

Before the development of FIT2D, generating the DT NS2D RANS solver input files 
was time consuming and tedious. A requirement was established for the development 
of a computer program that could quickly generate the flow domain and RANS input 
files for any 2-D hydrofoil that could conceivably be tested in the MIT-MHL. This 


stated need was translated into the following set of functional requirements or program 
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Figure 3-3: Typical D'TNS2D zonal boundaries for a foil in a tunnel 


objectives: 


ПЕ 


. Support basic inputs such as angle of attack, domain extents, and locations of 


Te 


Run from within a simple user interface. 


walls if any. 


. Allow the user to specify grid spacing in key areas such as the leading edge and 


trailing edge and inside the boundary layer. 


. Evaluate user spacing of elements inside the boundary layer and compare against 


flat plate calculations. Make à recommendation for any changes. 


. Read in z, y foil offsets. Spline the offsets and locate the leading edge. 
. Rotate the foil to user prescribed angle of attack. 
. Establish DTNS zone boundaries for a six zone “H” grid. 


. Allow for the zone boundary going downstream from the trailing edge to follow 


an established wake or a first guess wake from a vortex lattice calculation. 


. Generate all interior grid points for the mesh using isoparametric interpolation. 


. Write datafiles for the program TECPLOT containing a representation of the 


grid that can be viewed for correctness. 


Write datafiles for the INMESH grid smoothing program. 
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3.3.2 Structure of FIT2D 


FIT2D is written using the FORTRAN 77 code standard. It is executed interactively 
from a UNIX command window. The typical setup when running FIT2D is to have the 
data input file or control file displayed in an EMACS window(or any other editor), and 
a TECPLOT graphics window for viewing meshes. The three window configuration 
was chosen in favor of writing a dedicated X-Windows application. A UNIX command 
window, EMACS and TECPLOT provide an environment to which many people are 
already accustomed. After setting up the initial foil offset file(fname.foil) and control 


file(fname.ctrl), the user loops through the following basic procedure: 
1. Initialize FIT2D 
. Follow prompts to generate a grid. 
. View the grid in the graphics window. 
. If the grid is satisfactory, confirm permission to write INMESH files. 


сл >. CQ N 


. If the grid is unsatisfactory, make changes to fname.ctrl file and reinitialize 


BHD. 


Once satisfactory results have been achieved with FIT2D, the program INMESH is 
run to smooth the grid. This process is described in greater detail in section 3.6. 
Once the grid has been generated and smoothed, the geometry is ready for input to 
the RANS solver. A complete listing of the Program FIT2D is contained in Appendix 
À. 


3.3.3 Selecting Required Input 


Sample input files for FIT2D are contained in Appendix B. With the exception of 
input/output control, and boundary layer resolution changes, all FIT2D input is in 
batch form. The control file contains all of the data for describing the flow domain 
and how the mesh of elements will be distributed. The required input stems from the 
ideas presented in section 3.2.1. In the interest of making grid generation a somewhat 


more mechanical and less artistic process, it was decided early on to limit user input to 
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essential parameters for grid generation that the user must specify. 


Table 3.1: Batch Input Variables for FIT2D 


[ FOILGEO | frame foil geometry file name, located in same directory] 
[AOA | Angle of Attack in degrees 7) 
| Xpiv | distance z/c of pivot point from leading edge] 
| Ypiv [distance w/eofpivotpointfrommosetalline — — | 
| TUNPARAM | foil scale parameterc/tumndl wid —— | 
LUSL [upstream limit of flow domain in chord lengths from LE. — | 
[DSL —  |downstream limit of ow domain in chord lengths from TE. | 
[PHI forward rake distance of vertical zone lines from leading edge | 
[PAR [aft rake distance of vertical zone lines from trailmg edge | 
[NUS [number of elements in a-direction upstream  —  — | 
[NDS [mmber of elements in z-direction downstream | 
[Nys | mmber of elements in y-direction above and below fol | 
[NTOP | number of elements along chord on top surface of foil — | 
ГХВОТ | number of elements along chord on bottom surface of | 
|RESLE | Element width in chord lengths at leading edge | 
| RESMID | Element width in chord lengths at mid chord — — | 
[RESTE | Element width in chord lengths at trailing edge | 
[RESWALL | Element height in chord lengths at the walls | 
|RESBL | Element height in chord lengths on the foil surface | 
[PACK | fraction of NVS elements targeted within boundary layer | 
[Heg | chord based Reynolds number <  — | 
[ MESHFILE | INMESH input filename — č — | 
| RESTFILE | INMESH restart filename —— | 
[NUMIT —  [numberofINMESH ieraions — —  — — ^ —  — | 
| TOL | convergence tolerance limit for INMESH — — — | 


NRANSWK wake adaption flag: < O:straight line, 0: use VLM, > 0 use RANS data || 


| RANSWKPTS | z, y coordinates of the RANS wake line data 

































a few key parameters to guide the program. Within the control file, there are fifteen 


Additionally, 


there are some required administrative inputs. In the control file, each line of input is 


preceded by a descriptor line. The key inputs and their functions are summarized in 


table 3.1. Most of the parameters are intuitive. However, two parameters in particular 


need further explanation. The PHII and PHI2 rake angles are required in the mesh 


so that elements at the leading and trailing edge do not compress to a zero volume. 


An example of a problem leading edge is one that is circular with a large radius. If 


the zone boundary extended vertically from the leading edge, then the cells to the 


right of the boundary and above the leading edge would be skewed almost ninety 


degrees. In other words, cell volume is zero making it a singularity point. Referring 
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to Figures 3-1 and 3-3, one can see that the zonal boundaries are raked. The visible 
effect is that the cells retain a rectangular profile as they wrap around the leading 
edge. The upstream cells experience a minimal volume compression. But, the overall 
effect is to improve the consistency of the cell volumes. The raking capability was 


added at the trailing edge to accommodate highly curved anti-singing trailing edges. 
3.3.4 Splining and Element Distribution Methods 


All of the zonal boundary lines are splined parametrically in arc length using a cubic 
splining routine. The cubic coefficients are determined and stored for later use. Cubic 
splining in arc length was not required for any of the straight line boundaries. It was 
used anyway because it allows for future growth of the computer code to easily accept 
zonal boundary lines of arbitrary shape. An example of this is using a “C” shaped 
zone that wraps around the foil. 

The user specifies how many points are to be distributed along each boundary as 
well as the endpoint element dimension. Several schemes were explored to distribute 
the remaining points along the boundaries. The rate of change in cell dimensions 
should be relatively uniform along each boundary. Uniform spacing obviously is not 
a candidate. Cosine spacing was investigated, but it proved inadequate because the 
growth rate of cells near the areas of key interest was too fast. Ultimately, two 
spacing methods proved preferential for distributing elements. One spacing method 
is for horizontal zone lines. The other is for vertical zone lines. 

The horizontal lines define the tunnel walls, foil surface, upstream zone lines, and 
the wake. Along horizontal lines a parametric cubic distribution routine is employed. 
This routine was developed by Black[5] for distributing points on axisymmetric bodies 
and ducts. The program has been converted to a subroutine for inclusion in FIT2D. 
It has proven to be adequate for a variety of geometries. 

The spacing of elements on the vertical lines is somewhat more difficult. Vertical 


lines define how the elements are distributed normal to the foil and the tunnel walls. 
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There are two driving factors. One is to place an adequate number of very small 
elements close to the foil and walls to accurately capture the boundary layer profile. 
The other is to have much larger elements outside the boundary layer to minimize 
the total number. The user sets the spacing by specifying values for the variables: 
RESBL, RESWALL, NVS, and PACK. RESBL and RESWALL should be set using 
the y? method of section 2.7.4. Poor selection of RESBL and RESWALL (i.e. too 
large) will yield poor results and the boundary layer characteristics will be missed. 
The next trade-off is with NVS and PACK. It is always easy to increase number of 
points NVS, but the solution time suffers. PACK dictates how many of the total 
points NVS will be allocated to the boundary layer. 

Given that the steepest velocity gradients occur nearest to the surface, the fol- 
lowing spacing routine for vertical lines has been chosen. Set first element height 
at RESBL or RESWALL. Then geometrically grow the cells by a factor of 1.25 un- 
til the number of elements corresponding to PACK are used. Lastly, distribute the 


remainder of the cells using the cubic distribution routine. 
3.3.5 Defining Domain Boundaries 


Defining the length limits of the discretized domain was approached in several ways. 
For analyzing cases for foils in the tunnel, the actual physical limits of the tunnel 
were used to constrain the upstream and downstream limits of the gridded geometry. 
In the initial formulation of FIT2D, it was not foreseen that there would be a large 
need for analyzing unbounded flow around foils using the RANS code. However, as 
the research process progressed, more and more cases were analyzed in unbounded 


domains. 


The bounded case methodology. 


Appendix C shows the configuration of the MIT-MHL water tunnel test facility. 


Figure C-2 highlights the specifics of the test section where foils are mounted for 
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measurement. Once the scale of a subject foil is known with respect to the tunnel di- 
mensions, the wall locations and upstream and downstream limits can be calculated. 
These parameters are supplied to the input files as: USL, DSL, and TUNPARAM. 
FIT2D then places all boundaries at the proper location to replicate the tunnel ge- 
ometry. This methodology proved to be adequate. Results obtained agree well with 


MIT-MHL tunnel experimental results for lift, drag and wake profiles. 


The unbounded case methodology. 


For unbounded analysis the domain limits were extended. The boundary condition 
at the walls was changed from no-slip to tangent and the walls were moved far away 
from the foil. A convergence study was conducted using the vortex lattice subroutine 
to determine how far to locate the walls from the foil such that the flow around the 
foil could be considered unbounded. A distance of five chord lengths was established 
for the walls. At this distance, the difference between the lift coefficient calculated 
by VLM with images and without is approximately 0.1%. Pressure distributions on 
a foil surface were compared against results from the computer code PAN2D and the 
differences were negligible. The flow domain was extended upstream and downstream 
as well. Five chord lengths was judged suitable. Observation of pressure and velocity 
contours from the RANS outputs showed good uniformity at the inflow and outflow 


ends of the domain. 
3.3.6 Isoparametric Interpolation 


Once the zone boundaries have been identified and the cell spacing on the boundary 
is established, the initial location for the interior cell vertices must be established. 
This is a critical step in the grid generation process. The challenge is to locate all 
the interior points such that none of the quadrilaterals overlap and all the points 
are physically interior to the zone boundary. For a simple geometric shape such as 


a rectangle or a trapezoid the task is not overly demanding. Difficulties arise when 
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adjacent boundaries are curved or highly skewed. To overcome these difficulties, a 
methodology was devised wherein the arbitrary zone boundary geometry 1s paramet- 
rically mapped to a unit square. The interior points are determined using simple 
line geometry. Then, the interior points are mapped back out using the inverse map- 
ping function. This process is schematically shown in Figure 3-4. Numerically, this 
is accomplished using the subroutine INTERP which is a variation of isoparametric 


interpolation adapted from Bathe’s method for finite element formulation|4]. 
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Figure 3-4: Isoparametric Mapping Schematic of Interior Zonal Points 


3.4 FIT2D Output 


3.4.1 Tecplot Files 


Two ASCII plotting data files are generated by FIT2D. These files conform the data 
format required by the Amtec program TECPLOT version 7. 


The first file is rotate.plt. It contains the following: 
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e input geometry rotated to the desired angle of attack 
e an overlay of splined foil surface points used for foil surface zone boundaries 


e a horizontal reference line that passes through the z,y location of the foil pivot 


point. 


The purpose of this file is to validate the input foil geometry and to ensure that the 
pivot point and angle of attack have been specified properly. A sample of a figure 


created with rotate.plt is contained in Figure 3-5. 
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Figure 3-5: Display of rotate.plt data 


The second file is mesh.plt. This file contains the complete grid geometry for all of 
the zones. The purpose of this file is to allow the user to examine the grid computed 
by FIT2D to ensure that the desired spacing has been obtained and that the grid 
lines are smooth. Examples of plots generated with this file are presented in Figures 


3-1, 3-2, and 3-3. 
3.4.2 INMESH Input Files 


Once a satisfactory grid geometry has been obtained, the main output files from 


FIT2D are written. These files contain all the data required for running the computer 
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code INMESH(section 3.6). They are called the input and restart files. The input file 
is an ASCII data file which defines the zonal boundaries and how the cell points are 
distributed on the boundaries. The restart file is a binary formatted data file which 
contains all the interior zone point coordinates in addition to the zonal boundary 


points. 
3.5 Adjustment of Zone Boundaries in the Wake 


In section 2.4 the importance of aligning the grid zone boundaries aft of the foil to 
the convected wake was discussed. A special treatment was required to adapt the 
fluid mesh geometry to the expected wake line. Two methods were developed for 
wake adaption. The first is a vortex lattice method which finds the invicsid wake 
streamline that convects downstream from the trailing edge. A beneficial by-product 
of the vortex lattice solution is that a good approximation of the inviscid lift coefficient 
is obtained. The second method uses an initial RANS solution to extract a viscous 


wake streamline. 
3.5.1 Implementation of the ADAPT Subroutine 


A vortex lattice subroutine was implemented in the program which extracts the mean 
camber line of the subject foil in order to solve for the required circulation distribu- 
tion. Once the circulation distribution is known, the wake can be calculated. This is 
accomplished by extracting the field point velocity at the trailing edge and then incre- 
mentally marching in the direction of the field point velocity vector while extracting 
the field point velocity at each increment. This technique is visually demonstrated 
in Figure 3-6. The numerical implementation of the vortex lattice method in the 


subroutine ADAPT is a direct extension from the theory presented in section 2.4. 
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Figure 3-6: Vortex Lattice Geometry with Extracted Wake Line 


3.5.2 Adding Image Vortices 


The ADAPT subroutine was initially constructed without image vortices which would 
account for the wall effect. This proved satisfactory for analysis of unbounded flow 
around foils where there was little or no separation. However, the wake from the 
unbounded VLM solution tended to continue to convect away from the true wake 
when walls were present. Initially, it was thought that a correction could be applied 
to the unbounded VLM which would straighten out the wake. Several attempts were 
undertaken to develop schemes which redirected the velocity vector based on wall 
proximity and inviscid lift coefficient. This was made to work well for a single foil 
at a limited range of angle of attack. But the method was not robust for universal 
application and would require perturbing the parameters with each application. The 
addition of vortex images was then explored. Following the development in section 
2.4.3, pairs of image vortex lines were added to the ADAPT routine. The final 
results demonstrating the improved wake tracking are shown in Figure 3-7. The 
inset panel clearly shows how the unbounded solution continues to convect downward 


downstream while the wake of the VLM with images eventually becomes parallel to 
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the tunnel walls. 
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Figure 3-7: Comparison of wakes for VLM(unbounded) to VLM with Images(bounded) 


3.5.3 Convergence of the Vortex Lattice Method 


Studying the convergence criteria of the VLM with images is a two-dimensional prob- 
lem. The goal of the convergence study is to determine a minimum number of com- 
ponents that must be modeled in the computer code in order to obtain a reasonable 
solution. First, one has to consider what is an adequate number of vortex panels to 
demonstate convergence for the lift coefficient. Next, one must consider how many 
pairs of images must be included such that the addition of more image pairs has no 
effect on the solution. The results of the convergence study are presented in Figures 
3-8 and 3-9. Based upon the trend of the two curves, the final configuration selected 
for the VLM solution employs 40 vortex panels and 16 sets of image pairs. This 
yields an error in inviscid lift coefficient of less than 0.04%. The percent error of lift 
coefficent is defined using the single precision numerically converged value for the lift 
coefficient with 640 vortex panels and 128 sets of image pairs as the reference. 


Structuring the grid zones aft of the foil using the VLM wake was useful in cases 
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Figure 3-8: Vortex Lattice Panel Convergence Curve 
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Figure 3-9: Vortex Lattice Image Pair Convergence Curve 


where tunnel walls are present and there is little or no separation occurring at the 
trailing edge. The VLM wake was not a good approximation of the true wake for 
foils that had extreme camber distributions at the trailing edge or that had separated 
flow. In these cases the viscous effects on the lift coefficient become more important 
and the inviscid solution has less validity. When the viscous effects start to become 
significant, the VLM wake overshoots the true wake in the same manner that the 


VLM without images overshoots in the bounded flow case. 


62 





3.5.4 Adapting Wake Zone Boundaries Using RANS Output 


Despite refinements to the vortex lattice method, a significant difference persisted 
between the wake predicted by the VLM wake and the wake profile obtained using 
the RANS solution. Therefore, a method was developed to extract the wake data 
from an initial RANS solution based upon the the vortex lattice wake estimate. This 
is accomplished using the following steps: 

1. Generate grid using the VLM computed wake. 

2. Run the RANS solver for a few thousand iterations. 
. Extract data for the wake centerline streamline using TECPLOT. 
. Condition the data with the computer code GETWAKE. 
. Update the fname.ctrl file with the wake data. 


. Generate a new grid using the RANS wake data. 


Ny C» Сл A фр 


. Run the RANS solver to convergence. 


The program GETWAKE is listed in Appendix D. It is a simple conditioner that 
extracts a representative truncated set of datapoints from the TECPLOT streamline 
data. There is a major drawback to this wake correction process. It is very time 
consuming because the RANS solver must be used iteratively. It should be restated 
that this is only necessary in cases where separation has occurred near the trailing 
edge or the camber distribution on the aft part of the foil is unusual or both. A 


comparison between the two wake adaption schemes is shown in Figure 3-10. 


3.6 Using a Poisson Solver to Refine Grid 


3.6.1 The Program INMESH 


The computer code INMESH is used to refine the grid geometry output from FIT2D 
and generate the geometry input files for DTNS2D[8]. INMESH is a robust grid 
generation and refinement tool that has been used in a variety of fluid mesh geometry 
applications. INMESH uses an elliptic solver to approximate the Poisson equation 


over a prescribed domain with specified boundary values. The calculated streamlines 
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Figure 3-10: Sample Comparison of VLM Wake and RANS Wake Zone Adaption 


and equipotential lines are used to formulated the grid. The forcing term in the 
Poisson equation is used to cluster the grid near surface boundaries and to ensure 
orthogonality. 

INMESH uses the Successive-Over-Relaxation( SOR) method to solve the system 
of equations. The discretized equations are solved by iteratively sweeping through the 
domain and shifting the values at the nodes until convergence is achieved. A superb 


analogy for this process was described by Black[6]: 

The best analogy to this solution process would be to consider the initial 
guess at a grid to be a bedspring system where each node has four springs 
attaching it to its neighbors. The spacing of the nodes on the boundaries 
is controlled by the user and clustering could be seen as variable strength 
springs. The SOR method releases each node and allows it to move to its 
'natural' position. By iteratively sweeping through the domain, the nodes 
will eventually stabilize to the final grid. 

INMESH is capable of commencing the grid generation process with only the 
boundary points specified. Then begins starts with a set of uniformly spaced interior 
points and starts the bed spring process. By employing the isoparametric spacing 


routine in FIT2D, a reasonable first guess at the location of all the interior points 


is obtained. These points are provided to INMESH in the binary restart file. A 
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significant benefit from doing this is that the INMESH code requires several thousand 


fewer iterations to achieve the same convergence tolerance. In general, a higher 


relaxation coefficient can also be used. Figure 3-11 demonstrates the results of grid 


refinement using the program INMESH. 
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Figure 3-11: Typical Grid Spacing Before and After Smoothing with INMESH 


Correction of INMESH output for use in DTNS2D 


The INMESH code is used for a variety of grid generation applications. It provides 


the input files for the DTNS family of RANS solvers. INMESH writes out DTNS 


input files that overlap the individual zones by one cell at each adjoining boundary 


which are used if higher order boundary conditions are desired. The DTNS2D variant 
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of the code does not support this overlapping geometry configuration. Therefore a 
small computer code is used to condition the INMESH output so that it is compatible 
with the DTNS2D code. The program is called PATCH and it is listed in Appendix 


E. PATCH is specifically tailored for a six zone ^H" grid scheme. 


3.7 The RANS Solver 


3.7.1 Overview of DTNS2D 


The flow solver used in this thesis is the David Taylor Navier-Stokes Two Dimensional 
computer code developed by Gorski[13]. The DTNS2D computer code has been 
validated with experimental results for flow around two dimensional lifting bodies by 
Nguyen[25]. 

The DTNS2D code solves the RANS equations for incompressible turbulent flow 
which are derived in section 2.6. The flow domain is discretized into a large num- 
ber of quadrilateral cells. The solution method is a cell centered finite difference 
method based on a finite volume derivation. The program is structured so that the 
flow domain can be broken into a number of separate blocks each having its own set 
of boundary conditions. This allows for modeling of flow around foils in tunnels or 
unbounded domains. Zone junctions support a number of possible boundary condi- 
tions such as: no-slip, tangent, continuity, pressure specified, or velocity specified. A 
typical zone scheme and associated boundary conditions was demonstrated in Figure 
3-3. 

The code utilizes a third-order upwind difference total-variation-diminishing (TVD) 
scheme applied to the convection terms and a second-order central difference scheme 
applied to the diffusion terms. A first order Runge-Kutta forward time marching 


scheme is used to determine the steady state solution. 
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Chapter 4 


Data Post Processing Methods 


4.1 Introduction 


A variety of tools were used for post processing of the КАМ output data. At a 
minimum, it was desired to obtain the lift and drag coefficients for a foil at a given 
Reynolds number and angle of attack. For the purposes of comparing the RANS 
solution with other foil analysis computer codes, the pressure distribution on the foil 
surface was also required. During the study, other questions arose regarding specifics 
of the flow characteristics within the boundary layer on the foil surface and in the 
wake. For studying foils in tunnels, it is beneficial to be able to extract tunnel wall 


boundary layer profiles. 
4.2 UNNS2D 


The output from D'TNS2D contains all the flow characteristics at the center of each 
cell. However, the data is not in a readily viewable or workable format. A data 
conversion program UNNS2D, written by Scott Black of the MIT-MHL, was used 
to convert the data from raw output to a format compatible with TECPLOT v7.0. 
UNNS2D is a simple program and it could be easily adapted to convert data to any 
popular graphics software data format. The FORTRAN 77 source code for UNNS2D 
is located on the MIT-MHL software archive server. Once the data is loaded into 


TECPLOT, it can be displayed and manipulated in a variety of ways to view field 
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pressure and velocity contours or to draw streamlines. The UNNS2D output file is 
called dtns2d.dat. It contains data in columns by zone for each cell centered coordinate 


in the following order:z, y, u, v and P. 
4.3 Computing Lift and Drag 


Lift and drag are computed using the program FLD2D (Foil Lift and Drag Two- 
Dimensional). It is an adaptation of a program that was previously used in the MIT- 
MHL to calculate the lift and drag for the duct and body components on ducted 
axisymmetric bodies. FLD2D extracts the cell centered flow properties from the 
DTNS2D output files. The lift and drag are numerically integrated using the methods 
of section 2.3.3. The primary outputs from FLD2D are the shear component of the 
drag coefficient and the drag and lift coefficient obtained by integration of the surface 
normal pressure. The shear component of the lift is neglected because it is a very 
small quantity when compared to the lift obtained by pressure. À plot file called 
cp.dat contains data for the pressure distribution (Cp) on the foil and numerical 
trail for the integrated quantities. The FORTRAN 77 source code for FLD2D is also 


located on the MIT-MHL software archive server. 
4.4 Bounding Box Calculations 


One focal point of this thesis is to evaluate differences between computational results 
and experimental measurements. In the MIT-MHL Water Tunnel the CONTOUR 
bounding box program is used to calculate the lift and drag of 2-D foils from measured 
velocities around a bounding contour. The contour integration techniques used by 
the program CONTOUR are presented in section 2.5.1. A bounding box can be 
extracted from the dtns2d.dat file using TECPLOT and subsequently processed by 
the CONTOUR program. This is useful for a variety of comparisons. For instance, a 


side by side comparison of contour normal velocities can be made between identical 
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boxes used in an experiment and an equivalent geometry RANS solution. Bad data 
points in experimental measurements can be identified. Or if good experimental 
data is available, poor modeling in the RANS code or grid scheme can be identified. 
Additionally, bounding box calculations when applied to the RANS output data, 
provide a validation of the surface pressure integration routine(FLD2D) for lift and 


drag. 
4.5 Extracting Velocity Profiles in the Boundary Layer 


In the early stages of this thesis, research focused on the validation of previously 
obtained lift and drag characteristics of foils. It became apparent that the details of 
the flow inside the boundary layer near the foil and tunnel walls were also of interest. 
Using the LDV apparatus, it is possible to take velocity data at a sufficient number of 
points inside the boundary layer to get an accurate representation of the profile. The 
same data can be extracted easily from the RANS output using the TECPLOT data 
extraction feature. Once the profiles are extracted, characteristics of the boundary 


layers can be quantified using the methods of section 2.7. 
4.5.1 Post Calculation of y* 


A validation check of the grid adequacy in the boundary layer region is obtained by 
extracting RANS velocity and position data from the first several cells above the foil 
surface. A reverse calculation of the y* value at these three locations is performed. 
The y* values for cell spacing are compared with the guidance in section 2.7.4. The 
computer code YPLUS is the small program used for this task. The FORTRAN 77 
source code for YPLUS is located on the MIT-MHL software archive server. 

The results for lift and drag are highly dependent upon the spacing of the first 
several cells above a no-slip surface. To demonstrate this a foil was analyzed using 
the FIT2D/DTNS2D analysis tool. For each case, the spacing of the first several 


cells above the foil surface was varied. The maximum y* used for the first cell was 
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approximately twenty. The minimum y* value used was less than one. All other 
grid parameters were held constant. Figure 4-1 shows plots of the pressure coefficient 


for each case. The cases corresponding to larger initial y* values yield higher lift 
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Figure 4-1: Pressure Distribution Dependence on y* Parameter 


coefficients. The pressure distribution for the final case, with y* set at 0.8, compares 


well with other numerical results. 
4.6 Extracting Velocity Profiles in the Wake 


It was suspected that the turbulence models used in DTNS2D were inadequate. Data 
were extracted from the RANS solution in vertical cuts downstream of the trailing 
edge of the foil. These data were compared to identical wake cuts measured experi- 
mentally. Once again, TECPLOT was useful for extracting this data. These results 


are presented in Case Study I for the HRA foil. 
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Chapter 5 


Results 


In this chapter, two case studies are presented. The first is an analysis of a proposed 
advanced foil section at several angles of attack and a single Reynolds number. The 
flow around the foil is computed using DTNS2D in unbounded and bounded flow 
domains. Unbounded results are compared with the 2-D lifting surface analysis code 
PAN2D-BL. Bounded results are compared with experimental measurements. Differ- 
ences and similarities are highlighted. A simple correction scheme is presented which 
extrapolates bounded measurements for lift to unbounded values. 

In the second case study, a foil with a highly cambered trailing edge is analyzed at 
a single angle of attack and Reynolds number. The foil is analyzed using DTNS2D 
in unbounded and bounded flow domains. Unbounded and bounded results are com- 
pared with other computer codes and experimental measurements. This foil has 
proven to be very difficult to analyze with current numerical methods. The results 
of various computer codes are different enough that further detailed experimental 
study is warranted. The results in this case study are used to indicate regions of the 
flow around the foil where additional experimental data should be obtained. This 
case study demonstrates how RANS analysis can be useful for developing a stream- 
lined experimental test plan which will maximize the useful data obtainable while 
minimizing the costly resources and time associated with experimental study. 


These case study results are based upon DTNS2D solutions obtained from nu- 
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merical grids that demonstrated convergence with optimum y* spacing. The specific 


parameters used for the generation of these grids are presented in Appendix B. 


5.1 Validity of the Results 


5.1.1 Sources of Error 


There are several potential sources of error in the DTNS2D RANS results. There is 
no provision in the computer code for modeling the laminar regime of a boundary 
layer. DTNS2D uses a boundary layer model that starts out fully turbulent from 
the leading edge of the foil or the beginning of a wall. Along the walls this error 
is probably not significant. On the foil, this effect may or may not be significant. 
The real transition point from laminar to turbulent flow is different on the pressure 
and suction sides of the foil. The errors resulting from the omission of modeling 
the laminar boundary layer and transition are difficult to quantify. Intuitively, the 
magnitude of this error is most likely dependent upon the Reynolds number, the 
camber and thickness distributions and upon the angle of attack. 

Another phenomenon which may occur is that the turbulent boundary layer may 
revert back to a laminar flow after transition has occurred in certain regimes of 
Reynolds number. This effect is nearly impossible to quantify due to the stochastic 
nature of the process. But, nonetheless, it is another aspect of the real fluid flow that 


is omitted in the RANS solution which may have some impact on the overall results. 


5.1.2 General Comments Regarding Comparison with Tunnel Experi- 
ments 


In experimental tests, boundary layer transition on the foil is usually forced to occur 
on the foil at a specific location by the introduction of a turbulent stimulating surface 
irregularity such as bumps, rivets or a trip wire. At the location of the turbulence 
stimulators turbulent transition will occur with a fair degree of certainty. Forcing 


transition allows a more equitable comparison between experimental measurements 
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computer codes such as PAN2D-BL and XFOIL. By tripping the real flow in the 
experiment and specifying this same location in the computer code, one can consider 
the flows to be identical. 

Although the foil section mounted in a tunnel is two-dimensional, the flow within 
the tunnel still retains some three-dimensional effects. In à pure two-dimensional 
numerical flow solver, the side walls are completely omitted. In the real tunnel, the 
side walls tend to have uneven boundary layer growth on the areas above and below 
the foil. The magnitude of these effects are just now receiving attention and are being 
evaluated through related research in the MIT-MHL[19]. 

The RANS computer code models absolute steady state conditions that are impos- 
sible to achieve in a water tunnel. In the tunnel there are minor variations in inflow 
velocity. In the real flow unsteady vortex shedding can occur at the trailing edge. 
Very thin foils operating at high lift coefficients may twist and flex under the loading. 
The tunnel is built to standard engineering tolerances. The measured angle of attack 
in the tunnel can only be set to within several hundredths of a degree. While in the 
gridded model the geometry is precisely specified. All these issues combine and add 


to the resulting uncertainty of the experimental results. 
5.2 Case Study I: The HRA Foil 


A proposed advanced foil design created by Hydrodynamics Research Associates( HRA) 
was selected for this case study. A significant amount of experimental data has been 
recorded for this foil in various configurations in the MIT-MHL Water Tunnel test fa- 
cility. Therefore, it is a natural candidate for use in validation of the FIT2D/DTNS2D 
analysis tool. 

The HRA foil section, presented in Figure 5-1, is intended for application 1n pro- 
peller blade design. This foil is innovative because of its reduced camber distribution 


near the leading edge for shock free entry. The moderate camber aft allows the foil 
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to sustain a broad uniform pressure distribution on the suction side enabling it to 
operate at a substantial lift coefficient while minimizing the potential for cavitation 


inception. The foil offsets are contained in Appendix F.1. 
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Figure 5-1: The HRA Foil Shape 


5.2.1 Validation Check of Grid Adequacy 


As a first check of the validity of the results obtained from DTNS2D, the y* grid 
spacing and the boundary layer profiles were analyzed. Within each set of bounded 
and unbounded runs, the only parameter varied was the angle of attack. Therefore, 
only two different FIT2D control files were necessary: one for the bounded runs and 
one for the unbounded runs. The grids were identical except for the small variations 
associated with changing the angle of attack. Table 5.1 summarizes the post calcula- 
tion of the y* parameter for for the unbounded and bounded runs. The results shown 
are for the runs conducted at AO A = —0.28°. They are representative of all angles 
of attack that were analyzed. The first cell spacing falls within the target goal and 
those remaining conform to the distribution requirements specified in Section 2.7.4. 


The next check is to ensure that the boundary layer conforms to the Law of the 
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Table 5.1: Post Processing Check of y* for HRA Foil Grid 


[— — —Ls* Vale [vr Value] 
| Cell Number | Unbounded | Bounded | 
[Tagefri*|10 [ro] 
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Wall. The boundary layer profiles shown in Figure 5-2 are extracted at the mid chord 
of the suction side of the foil. These profiles are consistent over a large portion of the 
foil surface except near the leading and trailing edges. At the leading and trailing 
edges there is a significant amount of curvature and fairly steep pressure gradients. 
It is natural to expect deviation from the Spalding formula because it is based on flat 


plate theory. Based upon the y* spacing and B-L profile results, the first indication 
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Figure 5-2: Comparison of DTNS2D Computed B-L Profile with Spalding Formula for HRA Foil 


is that the results should be acceptable. The next step is to compare these solutions 


with other analysis methods. 


75 





5.2.2 Unbounded Flow Comparison 


In the unbounded case the lift and drag coefficients obtained from the RANS calcula- 
tions are compared with PAN2D-BL. A typical comparison of the foil surface pressure 
distribution is presented in Figure 5-3. This plot is prepared for a single case of the 
HRA foil at a 1? angle of attack. Comparisons at the other angles of attack eval- 
uated are similar. Throughout a range from —2° < a < +1° the FIT2D/DTNS2D 
analysis tool agrees with PAN2D-BL to within one percent of the lift coefficient. 
FIT2D/DTNS2D slightly over-predicts the total lift which is evident from the differ- 
ence in the pressure distribution near the trailing edge as shown in the inset panel 
of Figure 5-3. Values for drag coefficient computed by RANS and PAN2D-BL are 
comparable. Drag calculated from the RANS data is about 4% higher than the 
PAN2D-BL results throughout the angles of attack studied. A summary plot of these 


results is also presented later in Figures 5-6 and 5-7. 
5.2.3 Bounded Flow Comparison 


For the bounded case the lift and drag coefficients obtained from the RANS calcu- 
lations are compared with experimental results. Throughout a range from —2? «€ 
a < +1" the FIT2D/DTNS2D analysis tool agrees with the experimental results to 
within a few percent of the lift coefficient. FIT2D under-predicts lift compared to 
the experimental results. Some error may be caused by position measurement error 
in the water tunnel experiment. 

Values for drag coefficient computed by RANS and from experimental measure- 
ments are not as consistent as they are for the unbounded case. Drag calculated from 
the RANS data does not follow the same trend as the experimental measurements. At 
lower lift coefficients, such as AOA < —0.5°, RANS results agree well with the exper- 
imental results. There is a large jump in the RANS computed drag for higher angles 


of attack. This is inconsistent when compared with the trends of the unbounded 
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Figure 5-3: Comparison of Pressure Distribution on HRA Foil Using DTNS2D and PAN2D-BL(Re = 
25:90 104 — 1?) 


results and calculations. This may be caused by inadequate discretization of the flow 
domain or from a breakdown of the turbulence model under these flow conditions. A 


summary plot of these results is also presented later in Figures 5-6 and 5-7. 
9.2.4 Bounding Box Comparison 


Bounding box velocity contours extracted from the RANS data are used as input 
for the CONTOUR program. The results from the contour integration of the RANS 
data along with the surface integration results are presented in Table 5.2. The contour 
integration results compare favorably with the surface pressure integration results for 
all of the unbounded data. This indicates that the RANS solution is numerically 
consistent throughout the domain. In the bounded cases, however, the CONTOUR 
and surface pressure methods do not agree. It is not necessarily clear which method 


is in error for the bounded cases. 
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Table 5.2: CONTOUR and FLD2D Results for Lift and Drag Characteristics of the HRA Foil 
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A sample contour obtained from the DTNS2D results and an MIT-MHL exper- 
iment are shown in Figure 5-4. The V-n component of the fluid velocity around 
the contours are plotted for each leg of the contour. The contour normal velocities 
are nearly identical for the upstream, downstream and top contour legs. The bottom 


contour leg data does not agree aft of the mid chord of the foil. 
5.2.5 Wake Profile Comparison 


To make an observation of the wake characteristics aft of the foil, vertical cuts of 
data were taken at successive locations downstream of the trailing edge. This data 
was taken at the same locations both experimentally and from the computational 
results. These cuts are presented in Figure 5-5. The first profile is taken 0.04 chord 
lengths downstream of the trailing edge. The experimental data and the RANS data 
are nearly coincident at this location. This implies that at points very close to the 
trailing edge the RANS solver is accurately capturing the true flow conditions in the 
vicinity of the trailing edge. The successive cuts still show good agreement. However, 
it can be observed that the RANS calculated wake is diffusing at a faster rate than 
the experimental wake. This may be due to ineffective turbulence model application 


by the RANS solver in the wake region. 
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Figure 5-4: Comparison of Normal Velocity Component Around Contour Box for DTNS2D and 
HRA Experimental Results(Re = 3 x 10%, AOA — —0.636?) 


5.2.6 Relating Bounded Measurements to Unbounded Characteristics 


The unbounded and bounded results for lift are summarized in Figure 5-6. It is 
desirable to develop a simple scheme to relate measurements taken in the tunnel 
to their unbounded values in the absence of walls. There are several options for 
implementing such a scheme. The key decision regarding which type of scheme to use 
depends on whether you believe the experimental results or the numerical results to 
be more accurate. Rather than debate the merits and faults of the these results, it 
is assumed here, for illustrative purposes, that the experimental results are accurate 


and have been obtained with a high level of confidence. 
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Figure 5-5: Comparison of Wake Profiles for DTNS2D and HRA Experimental Results(Re = 3 x 
10°,AOA = —0.28°) 


Given this assumption, a method is required to correct the results obtained from 
the unbounded computer code estimates to unbounded true results. The method 
formulated here neglects side wall boundary layer effects. It is assumed the RANS 
solver is fairly accurate at capturing the magnitude of the relative effects between the 
bounded and unbounded cases. Referring back to Figure 5-6, each of the lift curves 
is linear in the regime studied. The unbounded results for DTINS2D and PAN2D-BL 
are close enough that the results can be considered equal. A least squares linear 
regression for the lift curves yields the results for lift slope and intercept presented in 


Table 5.3. Each lift curve can be represented with the following equation: 
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Figure 5-6: Summary of Results for Lift Coefficient vs. AOA for the HRA Foil 


Table 5.3: Lift Slope and Intercept Data for HRA Foil Calculations 


MHL Expt. 0.137 C/AOA |] 0.267 CLOAOA = 0° | 
DTNS2D Bounded 0.247 | 
| 





PAN2D-BL/DTNSZD(unb) 


CL =a x AOA +b, (5.1) 


where a and b are the respective slopes and intercepts. The difference in the slopes 
and intercepts of the bounded DTNS2D calculated curve and the experimentally 
measured curve represent the relative error between the estimated and actual result. 


Slope and intercept correction factors are formulated by: 


Qcf — Gezpt — ADTNS2D(bounded) 
бо x Dexpt nn DDTNS2D(bounded) (5.2) 
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Figure 5-7: Summary of Results for Drag Coefficient vs. AOA for the HRA Foil 


Once side wall boundary layer effects have been quantified, additional terms could be 
placed on the right hand side in equation 5.2. The correction factors are applied to 


the unbounded numerically derived lift curve in the following manner: 
(Ca) ee nd (Gund чя ас) X AOA m (биль Zu ber): (5.3) 


Equations 5.2 and 5.3 are applied using the values presented in Table 5.3. Figure 5-8 
shows the corrected lift curve compared with the previously obtained results. 

What is learned from this case study is that a RANS solver can be used as a liaison 
between bounded and unbounded flows. The FIT2D/DTNS2D analysis tool agrees 
well with current state of the art unbounded foil analysis tools such as PAN2D-BL. 
In bounded flow the FIT2D/DTNS2D analysis tool agrees well with experimental 
measurements. Therefore, to conduct a complete analysis for a given lifting flow, 


where comparison with experimental results are desired, the following methodology 
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Figure 5-8: Lift Curve With Correction Factors Applied for HRA Foil 


Is proposed: 


1. Based on the fact that FIT2D/DTNS2D and PAN2D-BL agree well for un- 
bounded flow, use PAN2D-BL to evaluate the foil in the range of angle of attack 
and Reynolds numbers that will be used in the experimental study. 


2. Select a few angles of attack and Reynolds numbers to evaluate using the FIT2D/DTNS2D 
tool for the bounded case. 


3. Obtain experimental results and compare the results to the bounded FIT2D/DTNS2D 
results. 


4. Using equations 5.2 and 5.3, compute the correction factors for the unbounded 
predictions and calculate the corrected lift curve. 

The above methodology should make efficient use of computational assets and result 

in a consistent comparison method between experimental and numerical results. The 

corrected lift curve provides an estimate of the error between the actual measured 

lift and drag of a foil and the numerically predicted unbounded lift and drag. The 


long term goal of the correction factor scheme is to provide motivation to improve 
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measurement and numerical prediction techniques until the correction factors are 
reduced to zero. Once the correction factors are reduced to zero then it can be 
concluded that the numerical computer codes have captured all the essential physics 


of the problem. 
5.3 Case Study II: Foil With A Cupped Trailing Edge 


The development of the “Cupped Foil” is interesting. The parent shape of the cupped 
foil is derived from a NACA 0016 thickness distribution foil with an a = 0.8 mean 
line and a maximum camber of 2.55%. The parent foil was given the name “B1” 
and is representative of a typical destroyer propeller blade section at r/R = 0.7. 
The trailing edge was beveled using the standard U.S. Navy formula for anti-singing 
trailing edges. 

Bloch suggested that it was possible to design a foil with a substantially higher 
lift to drag ratio by shifting the maximum camber distribution aft and cupping the 
trailing edge downwards|[7]. This conclusion was based on a series of numerical results 
from the computer code XFOIL. Subsequently, Jorde[16] conducted experiments to 
compare the B1 parent foil to the B1 modified with a cupped trailing edge. The 
modified B1 foil section with the cupped trailing edge is presented in Figure 5-9. The 
foil offsets are contained in Appendix F.2. 

Although the cupping at the trailing edge is not extreme to the human eye, it is 
proving to be a substantial challenge to the various array of computer codes available 
for numerical analysis. Currently there is no agreement among the computer codes 
that have been used to analyze the foil. Therefore, it is an interesting foil to use for 
this case study as a validation of the FIT2D/DTNS analysis tool. The first goal in 
this case study is to identify differences between the FIT2D/DTNS2D tool and other 
computer methods. This information will be used to provide guidance for future 


experimental research involving this foil as to where 1t would be productive to gather 
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detailed data. Once a suitable database of experimental data is obtained reflecting 
the actual conditions of the flow around the foil, then changes can be incorporated in 
the computer codes to make their methods more accurate. The difficulties associated 
with the unique geometry of the cupped B1 foil provide a valuable opportunity to 


improve the robustness of state of the art computational tools. 


0.5 0.75 1 
x/Chord 


Figure 5-9: 'The B-1 Foil With Cupped Trailing Edge Modification 


5.3.1 Validation Check of Grid Adequacy 


As with the HRA foil, the first check of the validity of the results obtained from 
DTNS2D is to check grid spacing and boundary layer profiles. The cupped foil was 
studied at one angle of attack and Reynolds number for unbounded and bounded 
calculations. T'wo different control files were necessary: one for the bounded run 
and one for the unbounded run. Table 5.4 summarizes the post calculation of the 
y* parameter for for the unbounded and bounded runs. The results shown are for 
the runs conducted at AOA = 0.5? and Re = 3 x 109. The first cell spacing falls 
within the target goal and those remaining conform to the distribution requirements 


specified in section 2.7.4. 
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Table 5.4: Post Processing Check of y* for B-1 Cupped Foil Grid 


[TF Value Гу Value 
Y (actual) [08 [08 — 
p mms us 
]a3 — E 
























The next check is to ensure that the boundary layer conforms to the Law of the 
Wall. The boundary layer profiles shown in Figure 5-10 are extracted at the mid chord 
of the suction side of the foil. These profiles are consistent over a large portion of the 
foil surface except near the leading and trailing edges. At the leading and trailing 
edges there is a significant amount of curvature and fairly steep pressure gradients 


so it is natural to expect deviation from the Spalding formula which is based on flat 


plate theory. 


Law of the Wall 
o Unbounded DTNS2D 
x Bounded DTNS2D 





0 + 
log, (У) 


Figure 5-10: Comparison of DTNS2D Computed B-L Profile with Spalding Formula for B1 Cupped 
Foil 
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5.3.2 Unbounded Flow 


Using the post processed DTNS2D output from the unbounded model, velocity and 
pressure contours can be viewed for the flow domain. Velocity contours near the 
trailing edge are presented in Figure 5-11. On this plot flow separation can be observed 
on the suction side of the trailing edge of the foil over the last three percent of 
the chord. This separation is de-cambering the flow at the trailing edge causing a 
significant reduction in lift coefficient from the inviscid value obtained from the VLM. 


The flow is also decelerating beneath the cupped region. 





0.5 0.6 0.7 0.8 0.9 


Figure 5-11: Velocity Contours Near Trailing Edge of Cupped Foil(unbounded flow) 


Pressure contours around the leading edge of the foil are presented in Figure 5-12. 
The concentric circles of pressure contour point to the leading edge stagnation point. 
An interesting observation of the pressure contours is that the value of the stagnation 
point pressure exceeds the pressure that is analytically possible. This could result 
from artificial compressibility introduced into the numerics of the RANS solver(an 
additional equation of state) to run incompressible solutions. Another possibility is 
simply numerical noise. The cell sizes near the stagnation point are very small. The 
calculations for those cells may be near the numerical precision limit of the CPU. 


This stagnation pressure error is observable in the data from all three of the RANS 
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solvers compared in this case study. 
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Figure 5-12: Pressure Contours Near Leading Edge of Cupped Foil(unbounded flow) 


For the analysis of the flow around the cupped foil FIT2D/DTNS2D results are 
compared with a RANS study by C.I. Yang of DTMB and results obtained by Kerwin 
using PAN2D-BL. A comparison of the calculated surface pressure coefficient for 
the three methods is presented in Figure 5-13. In general, the three methods yield 
comparable results. The largest differences occur in the leading and trailing edge 
regions. The difficulty in comparing these results relates to the difference in leading 
edge stagnation pressure. PAN2D-BL has a maximum pressure coefficient(Cp) of 
exactly 1.0, the maximum in the RANS solutions are approximately 1.1 for DT NS 
and 1.15 for C.I. Yang. Correcting for these differences may yield better alignment 
of the pressure coefficient curves. In any event, the differences in the trends at the 
leading and trailing edges are significant enough to warrant detailed experimental 
study in these regions in order to validate the predictions of the different computer 
codes. Comparison of lift and drag coefficient results are presented in Table 5.5. The 
inviscid lift coefficient obtained by the VLM is double the viscous value obtained by 


RANS and PAN2D-BL. The viscous effects on lift are significant with this foil. 
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Figure 5-13: Comparison of Pressure Distribution on Cupped Foil Obtained by Two Different RANS 
Solvers and PAN2D-BL 


Table 5.5: Comparison of Unbounded C; & Cp Calculations for B1 Cupped Foil (Re = 3 x 


106, AOA= 0.5%) 

[Method [0 [65 | 
[VIM 1140585 уа | 
| FIT2D/DTNS2D | 0.579 10.0102] 
[CI Yang RANS [0.505 |- | 
[PANZD-BL — [0472 |- | 
[XFOIL оов] 











5.3.3 Bounded Flow 


As is demonstrated for the unbounded flow, velocity contours near the trailing edge 
are presented in Figure 5-14. On this plot flow separation is again observed on the 
suction side of the trailing edge of the foil over the last three percent of the chord. 
In the bounded case, the separated flow region 1s wider and extends farther aft than 
in the unbounded case. This separation is de-cambering the flow at the trailing 
edge causing a significant reduction in lift coefficient from the value obtained from 


the VLM. As with the unbounded case, the flow is decelerating beneath the cupped 
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region. 





-0.1 
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Figure 5-14: Velocity Contours Near Trailing Edge of Cupped Foil(bounded flow) 


Pressure contours around the leading edge of the foil are presented in Figure 5-15. 
In similar fashion to the unbounded case the concentric circles of pressure contour 
point to the leading edge stagnation point. Án anomaly in the figure is the discontinu- 
ity of one of the pressure contour lines forward of the leading edge. Close inspection 
of the data file indicates the anomaly is caused by a TECPLOT interpolation er- 
ror across the zonal boundary. Additionally, as is observed in the unbounded case, 
the value of the stagnation point pressure exceeds the pressure that 1s analytically 
possible. 

For the analysis of the flow around the bounded cupped foil, FIT2D/DTNS2D 
results are compared with two different RANS studies by C.I. Yang of DTMB and 
Lafe Taylor of Mississippi State University. À comparison of the calculated surface 
pressure coefficient for the three methods is presented in Figure 5-16. For the bounded 
case, the FIT2D/DTNS2D results are in closer agreement with the results by Yang 
than in the unbounded case. 


Comparison of lift and drag coefficient results are presented in Table 5.6. There is 
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Figure 5-15: Pressure Contours Near Leading Edge of Cupped Foil(bounded flow) 


Table 5.6: Comparison of Bounded Cz & Cp Calculations for B1 Cupped Foil (Re = 3 x 10°, AOA = 


0.5°) 

[Method |o C» | 
[vom 12017 | n/a | 
[FIT2D/DINS2D_| 0.6413 | 0.0117 | 
[CI Yang RANS [0.605 [- | 
| Lafe Taylor RANS | 0.578 |- | 














a similar disparity between the inviscid results and the viscous results as is observed 
in the unbounded case. Once again, the viscous effects on lift are significant with this 
foil. 

As one would expect, the presence of the tunnel walls in close proximity to the 
foil affects the lift and drag characteristics of the foil. Figure 5-17 demonstrates the 
difference in the resulting pressure distribution on the foil. To keep the compari- 
son simple, only the FIT2D/DTNS2D results are shown. The tunnel walls cause a 
flow cambering effect to the fluid flow which increases the lift coefficient from the 


\ 


unbounded value. 
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Figure 5-16: Comparison of Pressure Distribution on Cupped Foil Obtained by Three Different 
RANS Solvers 


5.3.4 Tunnel Wall Boundary Layer Comparison 


In the bounded case, the presence of the Bl cupped foil alters the boundary layer 
growth on the upper and lower walls of the tunnel. This effect is greater on the upper 
wall. A plot of the velocity contours throughout the domain is presented in Figure 
9-18. On the upper wall, ahead of the foil, the boundary layer growth follows the 
typical trend for a flat plate. Above the foil the boundary layer thins. Aft of the foil 
the boundary layer grows rapidly. These changes in the boundary layer were initially 
observed in the RANS calculations and subsequently validated by MIT-MHL water 
tunnel LDV measurements. Boundary layer profiles for the upper tunnel wall are 
presented in Figure 5-19. The DTNS2D solution shows excellent agreement for the 
boundary layer characteristics ahead of the foıl. Agreement remains good above the 
foil. In the region above and just aft of the trailing edge, the DTNS2D boundary layer 
profile deviates markedly from the experimetal results. More detailed experimental 
measurements are required in this region to identify the exact location and conditions 


along the wall where the DTNS2D solution begins to deviate from the experimental 


92 












Foil 
Without Walls 
With Walls 










Res z3e6 
-0.8 AOA=0.5° 


0 0.25 0.75 1 


0.5 
x/Chord 


Figure 5-17: Comparison of Pressure Distribution on Cupped Foil for Bounded and Unbounded 
Flow 


measurements. 

Displacement thickness(ó*) and momentum thickness (0) integrals are computed 
at several locations along the upper wall. These results are pressented in Figure 5-20. 
From the leading edge to the upstream limit of the flow domain, FIT2D/DTNS2D 
agrees with the experimental results. Downstream of the leading edge, the RANS 
results differ significanlty from the experimental measurements. This may be a result 


of poor turbulence modeling in the regions of rapid pressure changes along the wall. 
5.3.5 Separated Flow 


Separation occurs near the trailing edge of the cupped foil at the angle of attack and 
Reynolds number used in this case study. The separated flow is clearly observable in 
Figures 5-11 and 5-14. Another interesting way to observe separated flow is through 
the use of streamlines. Streamlines are obtained in the TECPLOT software package 
by following velocity vectors in the domain from a specified starting point. Figure 5- 
21 shows streamlines obtained from the DTNS2D RANS solution. Separation occurs 


over the last 3% of the foil. Figure 5-22 obtained from the C.I. Yang results is 
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Figure 5-18: u Contours for Cupped Foil in Bounded Flow 


included for comparative purposes. The separation zone in Yang’s result is similar in 
length of the zone, but it is much wider across at the trailing edge. From the global 
perspective of the total flow, this difference is minimal but it is nonetheless present. 
The difference in the two RANS results indicates that detailed experimental data is 


required to validate the flow predictions at the trailing edge. 
5.3.6 Developing An Experimental Test Plan 


The cupped foil is a challenging foil to analyze numerically. There are many in- 
consistencies in the results among the different types of numerical solutions. More 
experimental data is required to validate the different aspects of the computer codes. 
Obtaining good experimental data is a very difficult task. Additionally, precious time 
and resources associated with experimental studies preclude obtaining data “every- 
where”. Here the RANS results can be applied to help identify the key areas that 
should be measured. Two benefits are realized. First, the possibility of taking data 
in areas of the flow that are not significant is reduced. Second, areas of interest that 
may be overlooked and not measured are reduced. 

The cupped foil case study yields the following guidance for the experimentalist 


to conduct a study of the B1 Cupped Foil: 
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Figure 5-19: Upper Tunnel Wall Boundary-Layer Profiles 


e The RANS solutions for the upper wall boundary-layer growth exhibit unusual 
characteristics. Boundary layer profiles should be obtained at several locations 


along the tunnel wall to confirm this. (see section 5.3.4) 


e Separation occurs at the trailing edge of the foil. The extent of separation is not 
consistent among the RANS solvers. Detailed velocity profiles are required over 


the last 3% of the suction side of the foil with additional wake cuts down stream. 


e The velocity contours on the pressure side of the foil near the “cup” are inter- 


esting. A moderate amount of data should be obtained in this region. 
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Figure 5-20: Upper Tunnel Wall Boundary-Layer Displacement Thickness and Momentum Thickness 


e At the stagnation point, the pressure values obtained by the RANS solvers is non- 
physical. The stagnation point should be located experimentally and pressure 


measurements in the vicinity of the stagnation point should be obtained. 


e The different numerical computer codes do not agree on the lift and drag char- 
acteristics of the foil. Lift and drag should be measured with statistical re- 
peatability using bounding box contours. The lift and drag coefficients should 
be calculated and the V - n component on each leg of the contour should be 


compared with the identical contours in the RANS domain. 


This provides the starting point for a detailed study. Based on the current numerical 
predictions, a lot of progress can be realized with the above test plan at a single angle 
of attack and Reynolds number. Once the above tasks are complete and the computer 
codes are validated, then the study can be expanded to a range of angles of attack 


and Reynolds numbers. 
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Figure 5-21: Streamlines Near Trailing Edge of Cupped Foil(DTNS2D Solution) 
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Figure 5-22: Streamlines Near Trailing Edge of the B1 Cupped Foil(C.I. Yang Solution) 
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Chapter 6 


Conclusions and Recommendations 


6.1 Conclusions 


The computer code FIT2D has been developed to rapidly generate the geometry for 
the fluid flow domain surrounding an arbitrary foil shape at a specified angle of attack 
in the MIT Marine Hydrodynamics Laboratory (MHL) water tunnel. This geometry 
is provided as input data for the RANS solver DTNS2D. A suite of software tools have 
been developed to provide post processing analysis to compare the RANS solution 
with other numerical techniques and experimental measurements. 

Through the use of the FIT2D/DTNS2D analysis tool, it has been shown that 
RANS solvers are quite useful for observing the details of the fluid foil around hydro- 
foils. The RANS results coupled with experimental data provide a good validation 
database for other numerical lifting analysis tools such as PAN2D-BL and XFOIL. 

The importance of the y* parameter should be emphasized again. The cells in the 
grid near a no-slip surface must be sized to follow the guidance in section 2.7.4 in 
order to capture the significant viscous effects within the boundary layer. Failing to 
follow this guidance will undoubtedly yield questionable RANS solutions. 

Adapting the fluid mesh zonal boundaries to the wake of the foil was somewhat 
successful. Alignment allowed better concentration of the finer discretization in the 
region of the steepest velocity changes. The wake adaption will offer significant ad- 


vantages if the turbulence modeling in the wake is improved. 
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FIT2D is robust and it can be applied to a variety of foil geometries. With minor 
modification is could be used for multicomponent systems such as yacht sail and mast 
assemblies. A version of FIT2D has been converted for use to generate grids to rep- 
resent axisymmetric ducted propulsor vechicles for use with the DTNS axissyemtric 
flow solver. 

Our knowledge of the details of two-dimensional fluid flow around hydrofoils is 
incomplete. The current lifting analysis codes do not yet achieve the results that 
minimize the performance risk in proceeding from design to production of advance 
foil section propellers without extensive model testing. It is apparant from the two 
case studies presented in Chapter 5 that improvements in propeller design computer 
codes can be achieved by first looking back to the details of two dimensional lifting 
flows. When closure is brought to the remaining details of the 2-D lifting analysis 
codes, the lessons learned can be incorporated back into the propeller design codes. 
Then advanced foil sections will be able to be incorporated in future propellers with- 
out having to undergo the monumental effort that was required in the Advanced 
Technology Demonstation for propellers at the Carderock Division, Naval Surface 


Warfare Center. 
6.2 Recommendations 


The following areas require further research in order to advance the current level of 


knowledge of two dimensional fluid flow around lifting surfaces: 


e Current turbulence models in computer codes do not work well in the wake 
region behind a foil. New models based on detailed experimental tests could do 
better at matching the rate of wake diffusion as the wake travels downstream. 
Improving the turbulence models will also improve the RANS solver’s predictions 


in the boundary layer near a no-slip surface. 
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Foils such as the B1 cupped foil are not easily modeled by current computer 
tools. More experimental study is required for these types of foils if their ad- 
vantageous lift/drag characteristics are to be exploited in the future. The ex- 
perimental results should be used to provide validation data for improved foil 


analysis computer codes. 


FIT2D should be modified so that it works better for generating grid geometry 
for unbounded flows. Current grid spacing algorithms cluster cells too finely at 


the outer regions of the flow domain. 


The DTNS2D RANS solver should be modified so that boundary layers start 


out laminar and then undergo transition to a turbulent boundary layer. 


The boundary layer growth on the side walls of a water tunnel is uneven. The 
magnitude of this effect is currently unknown. Detailed measurements are re- 
quired to quantify this effect. Additionally, it would be interesting to compare 
these measurements with a RANS solution for an equivalent three-dimensional 


domain. 


The issue of excessive calculated stagnation pressure needs to be addressed in 
DTNS2D. It is not known if the problem results from artificial incompressibility 
within the code or if the code is not adequately imposing the specified pressure 


boundary condition at the downstream extent of the domain. 
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Appendix A 


FIT2D Progam Listing 


PROGRAM FIT2D 


Сжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжж 


AQAANQAAAARAAN”A ARAANAAANAAANAANAANAANANRANDAKDNDNAANAANANAANAAARAnNDAAANARARARAARAADR 


Prepares all input files for program INMESH for 2D foils 
operating in the MHL tunnel at various angles of attack. 


Reads in foil geometry file, splines offsets 
splits upper and lower surfaces adjusts to 

angle of attack set in control file and prepares 
input and restart files for INMESH. 


2/97 JD 
Incorporated ADAPT subroutine to use vortex lattice lifting 
line to find wake dividing line to make grid wake adaptive. 


3/97 
Updated ADAPT subroutine to include image vortices due to walls 


3/97 
Updated input and adapt to accept grid line data for wake 
from rans output. 


WRITTEN BY: JOHN DANNECKER, MIT, 1996, Last Modified: 06 FEB 97 


жжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжж 


Program executive control file: 


* о э e o 89 9 9 © е ә еее ее еее ә % ә 9? re 


! NOTE: A CARRAIGE RETURN MUST BE PLACED AT THE END ! 
: OF ALL DATA FILES 


9 9 © o ө ә ее аа % 9 9 % е е ө 9 9 ее е е ее $9 * 9 9 э 9*9 * $9 * * $9 е е ә е «еее ә ее ә ® э 9? е © ө 


"fname.ctrl" 

The ".ctrl" file contains the governing information 

which specifies: 
foil file: name of .naca or .foil file 
AOA: foil angle of attack ref to free stream. 
Pivot point: point about which foil is inclined 
Scale Parameter: ratio of chord length to tunnel dim 
USL: upstream flow domain limit in chord lengths 
DSL: downstream flow domain limit in chord lengths 
PHI1: leading edge zone offset number of chord lengths 

rake zone forward at tunnel wall. 


PHIi:trailing edge zone offset number of chord lengths 
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€3 03 0203 0302 0320200 0C300303€00000000000600060400050»05505300000000000 


rake zone back at tunnel wall. 


Note: PHIi, PHI2 shall always be positive and these 


parameters shall not exceed 90% if USL and DSL 
respectively. 
Grid Resolution Parameters. 

NUS: number horizontal points upstream 


All params NDS: number horizontal points downstream 
shall be integer NVS: number vertical points above and 
and be n*8-1 below foil 

points. NTOP: number of points along top of foil. 


File Format: 


NBOT: number of points along bottom of foil 

RESLE: specifies smallest element size near 
the leading edge 

RESMID: midchord resolution 

RESTE: specifies smallest element size near 
the trailing edge 

RESWALL: specifies smallest vertical element 

size near a wall. 
RESBL: cell height within boundary layer 
REYNOLD: chord based reynold number 


Note: horizontal spacing along wall boundaries 
is is proportional to the spacing 
along the arc of zone boundaries 
ahead along and behind foil. 


fname.dat - desired inmesh input file name 
NUMIT: NUMBER OF INMESH ITERATIONS 
TOR: CONVERGENCE TOLERANCE FOR INMESH 
NRANSWK: -1 do staight wake line 

O do vortex lattice 

>1 use RANS data. 


FILE HEADER (limit 40 Chars) 

fname.naca or fname.foil 

AOA PIVOTX PIVOTY 

TUNPARAM Scale Parameter (Chord/Tunnel Section) 
USL DSL PHI1i PHI2 

NUS NDS NVS NTOP NBOT 

RESLE RESMID RESTE RESWALL RESBL, PACKFACTOR 
REYNOLD 

fname.dat 

NUMIT 

TOL 


жжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжж 
Foil Geometry Requirements: 


The leading edge of the foil is defined as the ZERO Station. 
After foil rotation, the Pivot point becomes the ORIGIN for 
all further calculations. 


"fname.naca" 


A "NACA" geometry file contains the Station, Camber 
and Thickness parameters for a foil starting at the 
leading edge marching aft to the trailing edge. All 
points shall be normalized values. If no leading edge 
radius is specified, leading edge curvature will be 
splined based on station data. Set Leading Edge 

Radius to 1.0 for a splined leading edge. 
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осо со со со со со о о со со о со о осо со со осо ag AanaAaAaAKnAKAANAAANAANAAANKAAAAAAQAAQAAAAQAAAAQAAARAQAAAAAARAAAQAQANQAAARQAARAARAAARAARAA 


File Format: 


FILE HEADER (limit 40 Chars) 

1 (integer, this flags program to NACA type geometry) 
## Chord Length 

## Number of Stations 

Thickness(tmax), Camber(ymax) 

Station(x/C) 1/2 Thickness(t/tmax) Camber(y/Ymax) 


Leading Edge Radius Multiplier 


"fname.foil" 


File Format: 


A "FOIL" geometry file contains the X,Y points that 
describe the surface of the foil. Foil data points 
shall be sorted so that the file marches from the 
trailing edge lower surface to the leading edge and 

then back to the trailing edge along the upper surface. 
If the two TRAILING edge points are NOT coincident, then 
fit2d.f will connect the two points with a straight line 
and specify an additional mesh zone extending from the 
blunt trailing edge downstream. 


FILE HEADER (limit 40 Chars) 

2 (integer, this flags program to FOIL type geometry) 
## Chord Length (for normalization) 

## Number of Data Points 

X Y Training Edge Lower Surface 


X Y Leading Edge 


X Y Trailing Edge Upper Surface 


Here is the Zoning Scheme required by INMESH and DTNS2D: 


Assemble the zones for input to transfinite interpolation. 


Here is how an indivdual zone is specified: 


Side 1 


Side 3 


en — — aum Gump ш mmm C— ср mum anum co A A A um md cum GENES NEED GIL C) amu m A 
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Сжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжж 


Єз СУ СУ 


IMPLICIT NONE 
INTEGER NPOINTS, NFTYPE,MIDPRES,MIDSUCT,NOSE 

INTEGER NUS, NDS, NVS, NTOP, NBOT, NTEMP,I,NI,NTEMP2,J,NTEMP3 
INTEGER NSUCTi1,NSUCT2,NPRES1,NPRES2,NPRESSURE,NSUCTION 
INTEGER NUMIT,NRESTART,NWHAT ,NLINE,NZONE,M,N 

INTEGER NTOPD,NVSD,NDSD,NBOTD 

INTEGER NRANSWK, K 

REAL AOA, PIVOTX, PIVOTY, TUNPARAM, USL, DSL, PHI1i, PHI2 

REAL CHORDL,FULLTHK,CAMBMX, RADLE 

REAL RESLE, RESMID, RESTE, RESWALL,RESBL,YUPWALL , YBOTWALL 
REAL TEMP2, TEMP1,TEMP3,TEMP4,PACKFACTOR,REYNOLD 

REAL TOL,E1,RF 

CHARACTER FIN*20, PLOTFLi*20, MESHFILE*20, RESTFILE*20 
CHARACTER FOILGEO*20, FOILDES*40 

CHARACTER CDUM*40, QUERY*10, JUNKTEXT*40 

CHARACTER PROMPT2*34 

REAL PI,TWOPI,ZERO,ONE,HALF 

PARAMETER (PI=3.14159,TWOPI=6.28318, NI=200, ZERO=0.0, ONE=1.0) 
PARAMETER (HALF=0.5) 


ARRAYS 


REAL A(3),B(3),C(NI),D(NI),E(NI), XI(NI),YI(NI),THCK(NI) 
REAL PRES1X(NI),PRES1iY(NI),PRES2X(NI),PRES2Y(NI) 

REAL SUCT1X(NI),SUCT1Y(NI),SUCT2X(NI) ,SUCT2Y(NI) 

REAL PRES1XS(NI),PRES1YS(NI) , PRES2XS(NI) , PRES2YS(NI) 

REAL SUCT1XS(NI) ,SUCT1YS(NI) ,SUCT2XS(NI) , SUCT2YS(NI) 

REAL PRESSUREX(NI), SUCTIONX(NI) 

REAL PRESSUREY(NI), SUCTIONY(NI) 

REAL UPLINEX (NI) , UPLINEY (NI) , DOWNX(NI) , DOWNY (NI) 

REAL DOWNXS(NI) , DOWNYS(NI) , UPLINEXS(NI) ,UPLINEYS (NI) 

REAL WALL1X(NI) ,WALL1Y (NI) , WALL2X (NI) , WALL2Y (NI) 

REAL WALLSX(NI) ,WALLSY(NI) , WALL6X(NI) , WALL6Y (NI) 

REAL WALL3X (NI) ,WALL3Y(NI) , WALL4X (NI) , WALL4Y (NI) 

REAL VERT1iX(NI),VERT1Y(NI),VERT3X (NI) , VERT3Y (NI) 

REAL VERT2X(NI),VERT2Y(NI),VERTAX(NI) , VERT4Y (NI) 

REAL VERTSX(NI),VERTSY(NI), VERT6X(NI) , VERT6Y (NI) 

REAL VERT7X(NI) ,VERT7Y(NI) , VERTSX(NI) , VERTSY(NI) 

REAL ZONEi1X(NI,NI),ZONE1Y(NI,NI),ZONE2X(NI,NI),ZONE2Y(NI,NI) 
REAL ZONE3X(NI,NI),ZONE3Y(NI,NI),ZONE4X(NI,NI),ZONE4Y(NI,NI) 
REAL ZONESX(NI,NI),ZONESY(NI,NI),ZONE6X(NI,NI),ZONE6Y(NI,NI) 
REAL DZONESX(NI,NI),DZONESY(NI,NI),DZONE6X(NI,NI) ,DZONE6Y(NI,NI) 
REAL DZONE3X(NI,NI),DZONES3Y(NI,NI),DZONE2X(NI,NI),DZONE2Y(NI,NI) 
REAL RANSWKPTS(2,NI) 

INTEGER IC(6,8,4) 


CORO OR ROR IO OIG I ROI IORI fe ae ae ae fe ae ae ae a de a ake ae ak fe ae aa ae 


C 
C 
C 
C 
C 
C 
C 
C 
C 


Open Control file fname.ctrl: default filemanme is foil2d.ctrl 


Seek Parameters describes above, use info for later output file. 


OPEN INPUT FILE 


PROMPT2 = *PROGRAM CONTROL FILE*//? (*//*foil2d.crt1*//*) =, 
WRITE (*,'(A,$)?) PROMPT2 

READ (*,?(A)?) FIN 

TSslEINCI:1).EQ.’ >) RIN = *foil2d.ctrl? 

WRITE (*,'(A)?) OPENING FILE: ^ // FIN 

OPEN (UNIT=4 ,FILE=FIN,STATUS=’OLD’) 
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C READ IN INPUT DATA 
READ(4, ' (A) ) CDUM 
READ(4,? (A) ) JUNKTEXT 
READ (4,’(A)’)FOILGEO 
READ(4,’(A)’) JUNKTEXT 
READ(4,*)A0A, PIVOTX, PIVOTY 
READ(4, ' (A) ) JUNKTEXT 


a 


Tunnel param for setting wall distance 


READ (4, *) TUNPARAM 

READ (4, › (А)? ) JUNKTEXT 

READ (4,*)USL,DSL, PHI1,PHI2 
READ(4,? (A) ) JUUKTEXT 
READ(4,*)NUS, NDS, NVS, NTOP, NBOT 
READ(4,’(A)’)JUNKTEXT 
READ(4,*)RESLE, RESMID, RESTE, RESWALL, RESBL, PACKFACTOR 
READ(4, (A) ) JUNKTEXT 
READ(4,*)REYNOLD 

READ (4,’ (A)? ) JUNKTEXT 
READ(4,' (A)  ) MESHFILE 
READ(4,? (A) ) JUNKTEXT 
READ(4,? (A)  ) RESTFILE 
READ(4,' (A) S JUNKTEXT 

READ(4,*) NUMIT 
READ(4,’(A)’)JUNKTEXT 

READ(4,*) TOL 
READ(4,’(A)’)JUNKTEXT 


C 
C Adding capability to read RANS wake points 
READ(4,*) NRANSWK 
с 
IF (NRANSWK.EQ.ZERO) THEN 
WRITE(*,*) ? xxx ? 
WRITE(*,*) '**NO RANS WAKE DATA, PROCEED WITH VLM*¥**’ 
WRITE(*,*) ? xxx ? 
ENDIF 
C 
с 
IF (NRANSWK.LT.ZERO) THEN 
WRITE(*,*) ? xxx ? 
WRITE(*,*) ?sxxxxxxx*STRAIGHT WAKE ЗЕГЕСТЕ)жжжжжжжжжжх › 
WRITE(*,*) ?'x**xNO VORTEX LATTICE SOLN AVAIL:sxsx**x*xxxx? 
WRITE(*,*) ? xxx ? 
ENDIF 
C 
c Test to see if there are points 
IF(NRANSWK .GT.ZERO) THEN 
WRITE(*,*) ? xxx ? 
WRITE(* ,*) ?**##*#**READING IN RANS WAKE рАТАжжжжжжжжж › 
WRITE(*,*)?*xx*x*x*xxNO VORTEX LATTICE SOLN AVAIL*#**##*4%? 
WRITE(*,*) ? xxx ? 
READ (4,*) ((RANSWKPTS (1,k) ,RANSWKPTS(2,k)) ,K=1,NRANSWK) 
€ write(*,*) ((RANSWKPTS(1,k),RANSWKPTS(2,k)) ,K=1,NRANSWK) 
ENDIF 
с 
С oe he fe fe fe oft fe fe afc fe fe fe fe eee fe fe he fe fe fe hee ee ake ale ade ade ade ade ade ade ade ade ade ale ale ade ale ade ad ale ad ad ali ad ale ode ade ade ole ade ЖЖ 
с жжжжжжжжжжжжжжжжжжжЕМ) OF INPUT FILE READS жжжжжжжжжжжжжж 
C Ae ade ole ade ole ole ode ale ade ae ade ole de oic okc ok ade ok okc ok ok ale okc ok okc ok ale ade ade oke ade ade ade ade ad ade fe ade ae ale a ade fe ade ode ali oke okc ok ok ole ole ae ae ale ale ade 
C 
C OPEN ECHO FILE TO VALIDATE INPUT FILE 
C 
с PROMPT2 = ’CONTROL CHECK FILE’//’ (’//’check.dat’//’) =? 
е WRITE (*,?'(A,$)?) PROMPT2 
E READ (*,?(A)?) FIN 
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ME FIN.: 1).EQ.? ?) FIN = *check.dat’ 
WRITE (*,’(A)’) ’OPENING FILE: ^" // FIN 
OPEN (UNIT=3 , FILE=FIN , STATUS=’ UNKNOWN’ ) 


WRITE(3,*)CDUM 

WRITE(3,*)FOILGEO 

WRITE(3,*)A0A, PIVOTX, PIVOTY 

WRITE(3,*)TUNPARAM 

WRITE(3,*)USL,DSL,PHI1,PHI2 

WRITE(3,*)NUS, NDS, NVS, NTOP, NBOT,RESLE, RESMID, RESTE, RESWALL 
WRITE(3,*)RESBL,PACKFACTOR , REYNOLD 

CLOSE(3) 

CLOSE(4) 

QC *k ok KK KH HK HK HF FF HF FF FF FF KK KK KK KK KK FF FF FF FF dd kk 


0000000000000 


OPEN FOIL DATA FILE FOR BASELINE GEOMETRY 


ОООО О 


WRITE (*,'(A)?) "OPENING GEOMETRY FILE: ^" // FOILGEO 
OPEN (UNIT=7 , FILE=FOILGEO,STATUS=’ OLD’ ) 
READ IN INPUT DATA 

READ(7,’(A)’)FOILDES 


С) 


NERHEE SPEGTRIES FOIL TYPE 
2 = X,Y GEOMETRY 
1 = NACA TYPE FORMAT 


n M 


су Су СУ СУ С) 


READ(7,*)NFTYPE 
READ(7,*)CHORDL 
READ (7,*)NPOINTS 


READ IN X,Y DATA 


су Су СУ 


IF (NFTYPE.EQ.2)THEN 
DO 5 I=1,NPOINTS 
READ(7 , *) XI (I),YI(I) 
5 CONTINUE 
ENDIF 


READ IN NACA FORMATTED DATA 


aaa 


IF (NFTYPE.EQ.1)THEN 
READ (7 , *) FULLTHK , CAMBMX 
DO 6 I=1,NPOINTS 
READ(7,*)XI(I),THCK(I),YI(I) 
6 CONTINUE 
READ (7 , *)RADLE 
ENDIF 
GEOSE(T) 


FAKE HF FH k ok oe o ok oe oe oe oe oe ok ck oe oe ok ok ok oe ok oe ok ok oe oe ce ok oe ok oe ok oe ok oe ok oe ce ke oe ok ok k k k k k k e e oe e o e e 


IF THE FOIL IS A NACA TYPE, IT NEEDS TO BE CONVERTED TO 
X AND Y DATA FOR FURTHER USE IN THE PROGRAM AND FEEDING 


IN TO INMESH. 


CHECK TO SEE IF IT IS A NACA FOIL, IF SO GO TO CONVERSION 
SUBROUTINE. 


оса ао ооо ооо 


IF (NFTYPE.EQ.1) THEN 
CALL NACACONV(NPOINTS,XI,YI,THCK,RADLE) 
ENDIF 
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OPEN ECHO FILE TO VALIDATE FOIL GEOMETRY FILE 


PROMPT2 = *GEO VALIDATION FILE*//* (*//*foilgeo.dat'//*) = ” 
WRITE (*,'(A,$)?) PROMPT2 

READ (*,'(A)?) FIN 

BENGEINOL:1)9EQ.' ")-FIN - "Фо реотдас” 

WRITE (*,'(A)?) 'OPENING FILE: ^ // FIN 

OPEN (UNIT=8,FILE=FIN ,STATUS=’ UNKNOWN’ ) 


WRITE(8, *) FOILDES 

WRITE(8,*)NFTYPE 

WRITE(*,*) ’input debugger’, CHORDL,NPOINTS 
99 FORMAT (F8.4,F8.4) 


ааазппаааааасаштао 


ak k k k k k k k k k k k k k k kk A OK KK KEK KKK KK KEE KKK k k k k k k k KEK KK KKK KK KK KK KKK KK 


€ DO 7 I=1,NPOINTS 

с WRITE(8,99) XI(I),YI(I) 
cat CONTINUE 

с WRITE(8,*)CDUM 

с CLOSE(8) 

C 

C 

с 


NORMALIZE FOIL DIMENSIONS IF REQUIRED 
IF (CHORDL.NE.ONE) THEN 
CALL NORMFOIL(CHORDL,NPOINTS,XI,YI) 
ENDIF 
с 


с 
C $ k k k k k k k k k k k k k k k ak ak k k k ak k ak FF FH FE EEE EEE k k k k 


С 
C SEND FOIL GEOMETRY TO BE ROTATED FOR ANGLE OF ATTACK 
С 
С 
с IF (AOA.NE.ZERO) THEN 
CALL ROTANGLE(NPOINTS,XI,YI,AOA,PIVOTX,PIVOTY) 
с ENDIF 
C 
C 
С 
с IF (A0A.EQ.ZERO) THEN 
с WRITE(*,*) 'AOA IS ZERO, ORIGIN DEFAULTS TO LEADING EDGE.’ 
с ENDIF 
C ok ale al ae le ae ale de ade ale ade ade ale ode ae ade ale ade ad ale ade ol ae ade ad ade ad al oke oke ok ok ok ade al ad ad oke oke oc okc oc oke k ak ade ae ae fe k k abe ae ade ale ae A k ak abe al ae a ad ale ae al ae k k k k 
C 
C PREPARE DATA FILE FOR TECPLOT TO VERIFY FOIL ROTATION 
С 
С CHECK IF USER WANTS CHECK FILE: 


WRITE(*,’(A,$)’) *CREATE A TECPLOT FILE TO VIEW ROTATION?<n>:? 
READ(*,?(A)?) QUERY 


debugger 
write(*,*) query 


о с СУС, 


IF ((QUERY.EQ.’Y’).OR.(QUERY.EQ.’y’)) THEN 
PROMBT2 = ’FOIL TECPLOT.FILE*//* (//tmotate.plt*//*) = ? 
WRITE (*,’(A,$)’) PROMPT2 
READ (*,’(A)’) FIN 
IF (FIN(1:1).EQ.’ ’) FIN = ’rotate.plt’ 
WRITE (*,’(A)’) "OPENING FILE: ^" // FIN 
OPEN (UNIT=9,FILE=FIN,STATUS=*UNKNOWN?) 
WRITE(9,*) "VARIABLES = X,Y’ 

WRITE(9,*) ’ZONE T="Input Geometry"’ 
DO 8 I=1,NPOINTS 
WRITE(9,*)XI(I),YI(I) 
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8 CONTINUE 

WRITE(9,*) "ZONE T="Reference Line"’ 

WRITE(9,*)-1.0,zero 

WRITE(9,*)zero,zero 

WRITE(9,*)ONE,ZERO 
€ CLOSE(9) 

ENDIF 

C ak ak ak ak ak ke ok ke kc ok ke okc ok ok okc ke ok ak ea ee ae ae ae ee жж oleo ok ok ee fee ee ee ea Ж 
SPLIT FOIL IN TO FOUR ARCS. SPLINE THE ARCS AND RETURN GRADED 
POINT ARRAYS TO DEFINE ARCS. THIS WILL DETERMINE THE MESH DENSITIES 
FOR EACH ZONAL AREA. 


IDENTIFY SPLITTING BOUNDARIES, FOILZONE returns the integer array 
position of the splitting points. 


СЗ О СУС СУ СУ СУ СЗ 


CALL FOILZONE(XI,NPOINTS,MIDPRES,MIDSUCT,NOSE) 


NOW TAKE THE FOIL OFFSETS AND SPLIT IN TO EIGHT ARRAYS FOR 
FEEDING IN TO THE FNSPLT SUBROUTINE. 4 each of X & Y points 


су С) Су СУ 


CALL ARCMAKER(XI,YI,1,MIDPRES,PRES1X,PRES1Y) 

CALL ARCMAKER(XI,YI,MIDPRES,NOSE,PRES2X,PRES2Y) 
CALL ARCMAKER(XI, YI,NOSE,MIDSUCT,SUCT1X,SUCT1Y) 
CALL ARCMAKER(XI,YI,MIDSUCT,NPOINTS,SUCT2X,SUCT2Y) 


ae ok k k ak k ak k ak k k ak k ak ak k ak k k ak k ak k k k k ok oe oce ok ok ok ok ok KE KKK KK KK KE KOK KK KK KK KKK KK KK KEK KKK KK EK KKK k k k 


SEND TO SPLINING ROUTINE 
Spline First Interval TE to midchord pressure side 


AO С) С) 


NTEMP is number of points want to get back from fnsplt 
NPRES1-2INT((float(NBOT)-1.0)/2.0) 


aa 


CALL FNSPLT(MIDPRES ,NPRES1,RESTE,RESMID ,PRES1X,PRES1Y,PRES1XS, 
+pres1YS) 


C 


debugger 


write(*,*) 'came back from arcmaker ok, nPRES1 is:',NPRES1 
write(*,*) 'Doing Aft Pressure Side’ 


write(79,*) ’ZONE’ 
do 4990 i=1,npresi 
write(79,*) presixs(i),presiys(i) 
4990 continue 


end debugger 
Now spline from midchord pressure side to LE 


0000300000000 с 


NPRES2 = NPRES1+2 

NTEMP2 is number of points going in to fnsplt 

NTEMP2 - NOSE - MIDPRES +1 

CALL FNSPLT(NTEMP2,NPRES2,RESMID,RESLE , PRES2X, PRES2Y,PRES2XS, 
*PRES2YS) 


o 


debugger 
write(*,*) ’Doing Forw Pressure Side’ 
write(*,*) ’number of points GOING IN TO fnsplt:’,ntemp2 
write(*,*) ’number of points from fnsplt:',npres2 
write(79,*) ’ZONE’ 
do 5000 i-1,npres2 
write(79,*) pres2xs(i),pres2ys(i) 


000000 OQ 
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c 5000 continue 

c end debugger 

с 

С жжжжжжжжжжжжжжжжжжпон do other half of foil. 

с Now spline from LE to midchord suction side 

C NTEMP is number of points want to get back from fnsplt 
C 


NSUCT1-INT((float(NTOP)-1.0)/2.0) 


o C 


NTEMP2 = MIDSUCT - NOSE+1 

CALL FNSPLT(NTEMP2,NSUCT1,RESLE,RESMID,SUCT1X,SUCT1Y,SUCTIXS, 
*SUCT1YS) 

debugger 


write(*,*) ’Doing Forw Suction Side’ 
write(*,*) 'number of points GOING IN TO fnsplt:',ntemp2 
Write(*,*) 'number of points from fnsplt:',nsucti 
write(79,*) 'ZONE' 
do 5010 i=1,nsucti 
write(79,*) SUCTi1xs(i),SUCT1ys(i) 
5010 continue 


Now spline from Midchord Suction Side to Trailing Edge 


0000000000000 ао 


NSUCT2 = NSUCT1+2 

NTEMP2 is number of points going in to fnsplt 

NTEMP2 = NPOINTS - MIDSUCT+1 

CALL FNSPLT(NTEMP2,NSUCT2,RESMID,RESTE,SUCT2X,SUCT2Y,SUCT2XS, 
ТЕШСІ275) 


о 


debugger 
write(*,*) "Doing AFT Pressure Side’ 
write(*,*) 'number of points GOING IN TO fnsplt:',ntemp2 
write(*,*) ’number of points from fnsplt:’,nsuct2 
write(79,*) ’ZONE’ 
do 5020 i=1,nsuct2 
write(79,*) SUCT2xs(i),SUCT2ys(i) 
5020 continue 
end debugger 


Sk He He He Hee oc oe oc o ok oko oe ok ook lalala da lefa dle alle dede eee kc okc oic okeokc ok okeoke okeoxke eoe oe eoe oe x 
Rejoin the pairs of arcs and make this the final foil geometry 


Ó001001000000000 о 


CALL ARCJOINER(NPRES1,NPRES2,NPRESSURE , PRES1XS , PRES2XS , PRES1YS 
+,PRES2YS , PRESSUREX , PRESSUREY ) 


CALL ARCJOINER(NSUCT1,NSUCT2,NSUCTION,SUCT1XS,SUCT2XS,SUCT1YS 
T2SUCT2YS, SUCTIONX,SUCTIONY) 
€ 
CR ka aaa aa a a ea ak ak a 
с Append the plot file with the splined data 
с 
с 
IF ((QUERY.EQ.’Y’).OR.(QUERY.EQ.’y’)) THEN 
WRITE(9,*) ’ZONE T="Pressure Side"' 
DO 9 I=1,NPRESSURE 
WRITE(9,*) PRESSUREX(I),PRESSUREY(I) 
9 CONTINUE 
write(9,*) *ZONE T="Suction Side"? 
DO 10 I=1,NSUCTION 
WRITE(9,*) SUCTIONX(I),SUCTIONY(I) 
10 CONTINUE 
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CLOSE(9) 


ENDIF 

C *k ok ade de ade ade ade He He He He ade ade ade ade ade ae ade ade ade ad al oje oie oc ocokc oc oc oc oc oc oic oc oe okcokc oc oc okc e e e ok oco oc oc ok oic ok okc okc okc oke oc oc okc oce oxcocokc oce oc oxeoe oe 
с Define upstream geometry of tunnel: 
E 
E First do horizontal line extending from leading 
E edge to forward end of tunnel. 

NTEMP3-10 
с write(*,*) ntemp3 


CALL TUNB(USL,SUCTIONX(1),SUCTIONY(1),UPLINEX,UPLINEY,NTEMP3) 
CALL FNSPLT(NTEMP3,NUS,RESLE,RESMID,UPLINEX,UPLINEY,UPLINEXS, 
+UPLINEYS) 


Qo ok ok EFT HF ade ale ale ade ade ale ale ale ale ade ade ale ade ale ale ade ale ade ade ale ale ale ale ale ade ade ale ale ade ade ale ade ade ale ale ade ale ale ale ade ade ale + 


Define downstream geometry of tunnel: 
First do horizontal line extending downstream from 
trailing edge to end of tunnel 


€ € е э €* е е е е е е е е Ф o €» * * 9 * 9 € ә е е е есе е е € 9 9 е 9* 9 9» @ @ 


Adapt the GRID boundaries to the wake based on VORTEX LATTICE 
SOLUTION. 


00000000 


IF(NRANSWK.EQ.ZERO) THEN 

CALL ADAPT(PRESSUREX ,PRESSUREY , SUCTIONX ,SUCTIONY ,DOWNX ,DOWNY , 
+ NBOT,NTOP,DSL,NTEMP3, TUNPARAM) 

ENDIF 


aa 
о 


if ranswk greater than O then there is rans data 


IF THE DATA IS BAD THEN NRANSWK WILL BE SET TO -1 AND 
A STRAIGHT WAKE WILL BE USED. 


абзаааса 


IF(NRANSWK.GT.ZERO) THEN 

CALL RANSWAKE(DSL,PRESSUREX (1) , PRESSUREY (1) ,RANSWKPTS , 
+ DOWNX , DOWNY , NTEMP3 , NRANSWK ) 

ENDIF 


If NRANSWK is set to -1 then it will do a straight wake. 


со с Су Су Су 


IF(NRANSWK.LT.ZERO) THEN 

write(*,*) rex? 

write(#*,*)’*****DEFAULTING TO STRAIGHT НАКЕжжжжжжжж) 
WRITE(*,*) жжж) 

CALL TUNB(DSL,PRESSUREX (1) , PRESSUREY (1) , DOWNX , DOWNY , NTEMP3 ) 
ENDIF 


c (NIN,NOUT,DS1,DS2,XI,YI,XO,YO) 


с Spline downstream line. 
CALL FNSPLT(NTEMP3 ,NDS,RESTE, RESMID , DOWNX , DOWNY , DOWNXS , DOWNYS ) 


aa 


debugger 


write(79,*) ’ZONE’ 

do 5030 i=1,nus 

write(79,*) uplinexs(i),uplineys(i) 
5030 continue 


0000000 
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write(79,*) 'ZONE' 

do 5040 i-1,nds 

write(79,*) downxs(i),downys(i) 
5040 continue 


ооооо 


CA ade ade de de ade ae ade ade ade ale ale ade ade ade ale ale ale ade ae ae ale ade ade ale ale ade ade ode e e ale ode ale ale ol od жж Жжжж Жжжж a ak ak ak e e e ale a ad od a e e al ad k ж жж ж 
с Locate the upper and lower walls of the tunnel. 
с 

CALL WALLFIND(TUNPARAM, YUPWALL , YBOTWALL ) 


C 
(Cade ade ale ale ade ade ole ade ode ale ale ale a ade ade ode a ade ade ade ade ae ade ae al af e e ad ad ad de ad abe al al ad al ale ale ae ale ae ae ae ale ade ade ale ade ale ale e e al ae abe ae ae ae ae ae ae ae ae ae a k 
с Assign XY Points to all sections of the walls 
€ 
C XBEG ,XEND, YOWALL,NPOINTS, XIN, YIN, XOUT, YOUT 
C>>>>>>>>>>> 
С forward upper 
с 
TEMP2-SUCTIONX(1)-USL 
TEMP1-SUCTIONX(1)-PHI1 
CALL TUNC(TEMP1,TEMP2,YUPWALL, NUS , UPLINEXS , UPLINEYS , WALL1X , WALL1Y) 
c 
с forward lower 
c 
CAEL TUNC (TEMP 1, TEMP2,YBOTWALL,NUS , UPLINEXS, UPLINEYS, 
+WALL2X , WALL2Y) 
C>>>>>>>>>>>>>>> 
с aft upper 
€ 
TEMP1=SUCTIONX(NTOP)+PHI2 
TEMP2=SUCTIONX(NTOP)+DSL 
CALL TUNC(TEMP1,TEMP2, YUPWALL, NDS, DOWNXS , DOWNYS , WALLSX , WALLSY) 
С 
с aft lower 
с 
CALL TUNC(TEMP1,TEMP2,YBOTWALL,NDS , DOWNXS , DOWNYS, 
+WALL6X , WALL6Y) 
C>>>>>>>>>>>>>>> 
с above foil 
€ 
TEMP1=SUCTIONX(1)-PHIi 
TEMP2=PRESSUREX(1)+PHI2 
CALL TUNC(TEMP1,TEMP2,YUPWALL,NTOP ,SUCTIONX,SUCTIONY, 
*WALL3X,WALL3Y) 
с 
с below foil 
€ 
CALL TUNC(TEMP2,TEMP1,YBOTWALL,NBOT ,PRESSUREX , PRESSUREY, 
+WALL4X , WALL4Y) 
с 


с 
PLIIIIITITITTTTETITITIIIITTITITIPIIIPIIIIIIIIIIIPIIIIIPPIPIEPIIIIIIIPIIIPTTII: 
с 
© It is very important to have the proper vertical spacing of 
€ cells within the boundary layer. Determine the y+ and compare 
с with the user RESBL and RESWALL. 

CALL FINDYPLUS(RESBL,RESWALL , REYNOLD , TUNPARAM, USL, DSL) 


с 
EA a k k k k k k k a k k a a k a e a e a a k a a k ak A ak a a ake a a e ak a e ak a ak ak ke ak a e ak ak ak ak ake ak akc akc ake ak ak akc ake ake a 


С Develop upstream vertical point spacing. Each line will consist 
с of three segments. Very fine spacing at the ends as per 

c RESWALL RESBL 

С 

c#*****%*DO ALL THE ABOVE LINES РТАЗТ*жжжжжжжж 

C 

с Upstream vertical line (above) 


ТІ 





TEMP 1=UPLINEXS (NUS) 

TEMP2=UPLINEYS (NUS) 

TEMP3=WALL1X (NUS) 

TEMP4=WALL1Y (NUS) 

CALL TUND(TEMP1,TEMP2, TEMP3, TEMP4,RESBL,RESWALL , 
*VERTi1X,VERTi1Y,NVS,PACKFACTOR) 


a a 


Leading Edge vertical (above) 


TEMP 1=UPLINEXS (1) 

TEMP2=UPLINEYS(1) 

TEMP3=WALL1X (1) 

TEMP4=WALL1Y (1) 

CALL TUND(TEMP1,TEMP2,TEMP3 , TEMP4,RESBL,RESWALL, 
+VERT3X, VERT3Y ,NVS , PACKFACTOR) 


Trailing Edge vertical (above) 


о 


TEMP1=downXS(1) 

TEMP2=downYS (1) 

TEMP3=WALLSX (1) 

TEMP4=WALLSY (1) 

CALL TUND(TEMP1,TEMP2, TEMP3,TEMP4,RESBL,RESWALL, 
*VERTSX,VERTSY,NVS,PACKFACTOR) 


Downstream Vertical Line (above) 


о 


TEMP1=downXS(nds) 
TEMP2=downYS(nds) 
TEMP3=WALLSX (nds) 
TEMP4=WALLSY (nds) 
CALL TUND(TEMP1,TEMP2,TEMP3 , TEMP4 ,RESBL,RESWALL, 
*VERT7X,VERT7Y,NVS,PACKFACTOR) 
с 
C s kk ale e kc o ok ke DO ALL THE BELOW LINES e ode od ox o ode ade al ale al ok ale x 
С 
€ Upstream Vertical Line (below) 
С 
TEMP1-UPLINEXS(NUS) 
TEMP2-UPLINEYS (NUS) 
TEMP3=WALL2X (NUS) 
TEMP4=WALL2Y (NUS) 
CALL TUND(TEMP1,TEMP2, TEMP3,TEMP4,RESBL,RESWALL, 
+VERT2X , VERT2Y , NVS , PACKFACTOR) 


Leading Edge vertical (below) 


aa an 


TEMP 1=UPLINEXS (1) 

TEMP2=UPLINEYS (1) 

TEMP3=WALL2X (1) 

TEMP4=WALL2Y(1) 

CALL TUND(TEMP1,TEMP2,TEMP3,TEMP4,RESBL,RESWALL, 
+VERT4X, VERT4Y , NVS , PACKFACTOR) 


Trailing Edge vertical (below) 


0000 


TEMP1=downXS(1) 

TEMP2=downYS(1) 

TEMP3=WALL6X(1) 

TEMP4=WALL6Y(1) 

CALL TUND(TEMP1,TEMP2,TEMP3,TEMP4,RESBL,RESWALL, 
+VERT6X , VERT6Y , NVS , PACKFACTOR) 


гү? 





Downstream Vertical Line (bleow) 


оосо 


TEMPi=downXS(nds) 
TEMP2=downYS(nds) 
TEMP3=WALL6X (nds) 
TEMP4=WALL6Y (nds) 
CALL TUND(TEMP1,TEMP2,TEMP3, TEMP4,RESBL,RESWALL, 
*VERT8X,VERT8Y,NVS,PACKFACTOR) 
С 
с 
CORR CRG ккк жж 
с 


с Assemble the zones AND DO isoparametric interpolation. 
с 
E 
с Here is the zone geometry: 
ЕР СЕН... ӛ8 “ее есе... 
с | | | | 
с | | | | 
с | zone 1 | zone2 | zone 3 | 
с ЕЕЕ 21а | | | 
с | [<======—==== Ж ан | 
с | | | | 
с | zone 4 | zone 5 |  гопе 6 | 
с | | | | 
c | | | | 
с s Lh e eos а |... 8 T D | 
с 
с 
С Here is how an indivdual zone is specified: 
e 
ce Rag a a ы ыы 
C | Side 4 | 
е | | 
€ | ZONE "n" | 
e Side 1 | | Side 3 
c | | 
с | | 
с | 514е 2 | 
c Е-Е OS i | 
€ 
€ 
ec 
€ /N 
c | Array reference to enter finite 
с | interpolation scheme. 
с J | 
с ЕС а. > 
€ T 
€ 
€ 
C ZONE ONE 
C 
CALL BORDERS(ZONE1X,ZONE1Y,VERT1X,VERT1Y,UPLINEXS ,UPLINEYS, 
+ VERT3X,VERT3Y , WALL1X, WALL1Y,NUS,NVS,NI) 
C 
C ZONE TWO 
C 
C 
C debugger 
С 
€ do 5050 i-1,nvs 
С write(81,*) vertb5x(i),vert5y(i) 
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c 5050 continue 


с 
с 


оса 00 су С) СУ С) ооо соо 


су Су СУ СУ 


С 


end debugger 
CALL BORDERS (ZONE2X , ZONE2Y, VERT3X , VERT3Y , SUCTIONX , SUCTIONY 
+,VERTSX, VERTSY,WALL3X,WALL3Y ,NTOP ,NVS , NI) 


ZONE THREE 


CALL BORDERS (ZONE3X, ZONE3Y, VERTSX, VERTSY, DOWNXS ,DOWNYS, 
+ VERT7X,VERT7Y,WALLSX,WALLSY,NDS ,NVS,NI) 

ZONE FOUR 

CALL BORDERS (ZONE4X , ZONE4Y , VERT2X , VERT2Y ,WALL2X , WALL2Y, 
+ VERT4X , VERT4Y , UPLINEXS , UPLINEYS ,NUS,NVS,NI) 


ZONE FIVE 


CALL BORDERS (ZONESX, ZONESY, VERT4X , VERT4Y , wall4X ,wall4yY 
+, VERT6X, VERT6Y , pressureX, pressureY , NBOT ,NVS,NI) 
ZONE SIX 


CALL BORDERS (ZONE6X , ZONE6Y, VERT6X , VERT6Y , WALL6X , WALLEY, 
+ VERT8X,VERT8Y,DOWNXS , DOWNYS ,NDS ,NVS,NI) 


Сжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжж 


C 
с 


6010 
6000 


6030 
6020 


6050 
6040 


Open Tecplot File for viewing Mesh. 


PROMPT2 = ’Mesh View TECPLOT FILE’//’ (’//’mesh.plt’//’) = ’ 
WRITE (*,’(K,$)’) PROMPT2,’:? 

READ (*,?(A)?) FIN 

IF (FIN(1:1).EQ.’ ’) FIN = ’mesh.plt’ 

WRITE (*,'(A)?) ?"OPENING FILE: ^ // FIN 

OPEN (UNIT=82,FILE=FIN,STATUS=’UNKNOWN’ ) 


topzones 
Write(82,*)'VARIABLES =X,Y? 
write(82,*)’ZONE T="ZONE i" I=’,Nus,’ J=’,NVS,’ F=POINT’ 
do 6000 j=1,nvs 
do 6010 i-1,nus 
write(82,*) zoneix(i,j),zoneiy(i,j) 
continue 
continue 
write (82,*) ’ZONE T="ZONE 2" I=’,Ntop,’ J=’,NVS,’ F=POINT’ 
do 6020 j=1,nvs 
do 6030 i=1,ntop 
write(82,*) zone2x(i,j),zone2y(i,j) 
continue 
continue 
write(82,*)’ZONE T="ZONE 3" I=’,Nds,’ J=’,NVS,’ F=POINT’ 
do 6040 j-1,nvs 
do 6050 i-1,nds 
write(82,*) zone3x(i,j),zone3y(i,j) 
continue 
continue 
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6070 
6060 


6090 
6080 


7010 
7000 


с 


bottom xones 


write(82,*)’ZONE T="ZONE 4" I=’,Nus,’ J=’,NVS,’ F=POINT’ 
do 6060 j=1,nvs 
do 6070 i=1,nus 
write(82,*) zone4x(i,j),zone4y(i,j) 
continue 
continue 
write(82,*)'ZONE T="ZONE 5" I=’,Nbot,’ J=’,NVS,’ F=POINT’ 
do 6080 j=1,nvs 
do 6090 i=1,nbot 
write(82,*) zone5x(i,j),zone5y(i,j) 
continue 
continue 
write (82,*)’ZONE T="ZONE 6" I=’,Nds,’ J=’,NVS,’ F=POINT’ 
do 7000 j=1,nvs 
do 7010 i=1,nds 
write(82,*) zone6x(i,j),zone6y(i,j) 
continue 
continue 
close(82) 


CA a a ale ade da ad ale KKK KK KK KK KKK KK KKK KK KKK KK KK KKK KKK KKK KK KEK KKK KKK KK KKK KK KKK KKK 


с 
с 


с 
с 


AAN aa 0000 


aa 


Generate INMESH input file: 


CHECK IF USER WANTS CHECK FILE: 


write(*,*) "RERRRRRRARRARARRARARRA RAR RARA RARARRARAARARARARO 
write(*,*) ? ? 

write(*,*) ’View the grid on TECPLOT, Check boundaries for’ 
write(*,*) ’wrapping. If the grid is not correct or if’ 
write(*,*) ’you want to make changes, you may abort gen-’ 
write(*,*) ’eration of the INMESH file.’ 

write(*,*) ’ т 

Write(*,*) "RERRARRARRARRARRARRARARARAARARA RARA RA RARA RARAR 
write(*,*) ? ? 


WRITE(*,’(A,$)’) ’DO YOU WANT TO WRITE THE INMESH FILE?<n>:? 
READ(*,?'(A)?) QUERY 


debugger 
write(*,*) query 


IF ((QUERY.EQ.'Y?).0R. (QUERY.EQ.^y?)) THEN 
WRITE (*,?'(A)?) "WRITING FILE: ^ // MESHFILE 
OPEN (UNIT=80, FILE=MESHFILE , STATUS=’ UNKNOWN’ ) 


FORMAT FIRST LINE 


E1 is the orthogonality at boundry 
E1=1.0 


Relaxation parameter 
RF=0.5 


Tell INMESH it is a restart 
NRESTART=2 


Unknown parameter 
NWHAT=1 
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98 


96 


c 
c 


WRITE FIRST LINE OF THE INMESH FILE 
format(f4.2,f12.8,17,f8.4,2i4) 
format(i3,7i3) 

format(i3,2i5,a20) 

if (numit.eq.1)then 

nrestart=1 

endif 

write(80,98) E1,TOL,NUMIT,RF,NRESTART , NWHAT 


сжжжжжжжжжжжжжжжжжжжжж 0 THE FIRST THREE 20МЕЅжжжжжжжжжжжжжж 


C 


сжжжжжжжжжжжугіъе out Boundaries for Zone i 


c 


с 
с 


с 


Be 
с 


write(80,96)1,nus,nvs,’ жжж7опе 1***’ 


Zone 1 Line 1 

NLINE=1 

Dera (1C(1,1,1),I=1,8)/1,5,0,0,1,1,0,0/ 
write(80,97) (IC(1,I,NLINE) , I=1,8) 


call zonelines(zone1x,zonely,nus,nvs,NLINE,80) 


Zone 1 Line 2 

NLINE-2 

ШАҒА СІС(1,1,2),І-1,8)/2,1,0,0,1,2,0,0/ 
write(80,97) (IC(1,I,NLINE),I=1,8) 


call zonelines(zoneix,zonely,nus,nvs,NLINE, 80) 


Zone 1 Line 3 

NLINE=3 

ШИК СТСС(І, І,3),І-1,8)/3,0,2,1,1,1,0,0/ 
write(80,97) (IC(1,I,NLINE),I=1,8) 


Zone 1 Line 4 

NLINE=4 

DATA (IC(1,1I,4) ,I=1,8)/4,6,0,0,1,4,0,0/ 
write(80,97) (IC(1,I,NLINE) ,I=1,8) 


call zonelines(zoneix,zoneiy,nus,nvs,NLINE,80) 


сжжжжжжжждопе 2 


с 
с 


ao 


Write out Boundaries for Zone 2 
write(S0,96)2,ntopti,nvs,’  ***Zone 2***? 


Zone 2 Line 1 

NLINE=1 
DABAM(1€C(2,1,1),1I121,9)/1,0,1,3,1,150,0/ 
write(80,97) (IC(2,I,NLINE) ,I=1,8) 


Zone 2 Line 2 

NLINE=2 

DASA (TC(2,I,2),1I121,8)/72,6,0,0,2,2,0,0/ 
write(80,97) (IC(2,I,NLINE),I=1,8) 


write(80,*)zoneix(nus-1,1),zoneiy(nus-1,1) 
call zonelines (zone2x,zone2y,ntop,nvs,NLINE,80) 


Zone 2 Line 3 
NLINE=3 
DATAM Со. т. 3), 1=1,8)/3,0,3,1,1,1,0,0/ 
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write(80,97) (IC(2,I,NLINE),I=1,8) 


с 
C Zone 2 Line 4 
NLINE=4 
DATA (IC(2,1,4),1=1,8)/4,6,0,0,2,4,0,0/ 
write (80,97) (IC(2,I,NLINE),I=1,8) 
с 
Wwrite(80,*)zoneix(nus-1,nvs),zoneiy(nus-1,nvs) 
call zonelines(zone2x,zone2y,ntop,nvs,NLINE,80) 
Cc 


с 
C 

C сжжжжжжжждопе 3 

с Write out Boundaries for Zone 3 
с 


write(80,96)3,nds+1,nvs,”?  *x**Zone Зжжж? 


С Zone 3 Line 1 
NLINE=1 
БЕ С(З,1,1),І-1,8)/1,0,2,3,1,1,0,0/ 
write(80,97) (IC(3,I,NLINE),I=1,8) 


с 
C 
с Zone 3 Line 2 
NLINE=2 
DATA (IC(3,1,2),1-1,8)/2,1,0,0,3,2,0,0/ 
write(80,97) (IC(3,I,NLINE),I=1,8) 
С 
write(80,*)zone2x(ntop-1,1),zone2y(ntop-1,1) 
call zonelines(zone3x,zone3y,nds,nvs,NLINE,80) 
C 
Е Zone 3 Line 3 
NLINE=3 
DENTS CICCS.I,39,121,8)/3,2,0,053493,0,0/ 
write(80,97) (IC(3,I,NLINE),I=1,8) 
с 
call zonelines(zone3x,zone3y,nds,nvs,NLINE,80) 
с 
€ Zone 3 Line 4 
NLINE=4 
DATA (1IC(3,1,4),I=1,8)/4,5,0,0,3,4,0,0/ 
write(80,97) (IC(3,I,NLINE),I=1,8) 
€ 


write(80,*)zone2x(ntop-1,nvs),zone2y(nus-1,nvs) 
call zonelines(zone3x,zone3y,nds,nvs,NLINE,80) 
бб rh hehehe ehe hh 
C 
Сжжжжжжжжжжжжж[)0 ТНЕ ЗЕСОМО ТНВЕЕ 20МЕЗжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжж 
С 


сжжжжжжжжжжж Write out. Boundaries for Zone 4 


с 
write(80,96)4,nus,nvs,’ ***Zone 4**x? 
с 
с Zone 4 Line 1 
NLINE=1 
DATA (IC(4,1,1),I=1,8)/1,5,0,0,4,1,0,0/ 
write(80,97) (IC(4,I,NLINE),I=1,8) 
с 
call zonelines(zone4x,zone4y,nus,nvs,NLINE,80) 
& 
C Zone 4 Line 2 


NLINE-2 
DATA (IC(4,1,2),1=1,8)/2,6,0,0,4,2,0,0/ 
write(80,97) (IC(4,I,NLINE),I=1,8) 
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call zonelines(zone4x,zone4y,nus,nvs,NLINE,80) 


c Zone 4 Line 3 
NLINE=3 
ЯНА (ІС(Ш І,3),І-1,8)/3,0,5,1,4,1,0,0/ 
Write(80,97) (IC(4,I,NLINE),I=1,8) 


с Zone 4 Line 4 
NLINE-4 
DATA (IC(4,1,4),1=1,8)/4,1,0,0,4,4,0,0/ 
write(80,97) (IC(4,I,NLINE),I=1,8) 
С 
call zonelines(zone4x,zone4y,nus,nvs,NLINE,80) 
Cc 
с 
сжжжжжжжж7опе 5 
е Write out Boundaries for Zone 5 
с 
write(80,96)5,nbot+1,nvs,?” жжждопе Бжжж) 


с Zone 5 Line 1 
NLINE-1 
DATA (IC(5,I,1),I=1,8)/1,0,4,3,4,1,0,0/ 
write(80,97) (IC(5,I,NLINE),I=1,8) 


€ 
С 
с Zone 5 Line 2 
NLINE-2 
А (ІС(Б.ІГ,2),121,87/2,6,0,0,5,2,0,0/ 
write(80,97) (IC(5,I,NLINE),I=1,8) 
с 
write(80,*) zone4x(nus-1,1),zone4y(nus-1,1) 
call zonelines(zone5x,zone5y,ntop,nvs,NLINE,80) 
€ 
с Zone 5 Line 3 
NLINE=3 
Паша нов T.3),I-1,8)/3,0,6,1,4,1,0,0/ 
write(80,97) (IC(5,I,NLINE),I=1,8) 
с 
с Zone 5 Line 4 
NLINE-4 
DATA (1C(5,1,4),I=1,8)/4,6,0,0,5,4,0,0/ 
write(80,97) (IC(5,I,NLINE),I=1,8) 
Е 
write(80,*) zone4x(nus-1,nvs),zone4y(nus-1,nvs) 
call zonelines(zone5x,zoneSy ,ntop,nvs,NLINE, 80) 
Ce 


с 

С 

C сжжжжжжжжопе 6 

Е Write out Boundaries for Zone 6 
с 


write(80,96)6,ndsti,nvs,’ жжж7опе бжжж? 


с Zone 6 Line 1 
NLINE-1 
DATA (IC(6. I,1),I=1,8)/1,0,5,3,4,1,0,0/ 
write(80,97) (IC(6,I,NLINE),I-1,8) 


зо 


С Zone 6 Line 2 
NLINE-2 
DARTE (6, 1,2),1=1,8)/2,6,0,0,6,2,0,0/ 
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Wwrite(80,97) (IC(6, I, NLINE),I-1,8) 


write(80,*)zone5x(nbot-1,1), zone5y(nbot-1,1) 
call zonelines(zone6x,zone6y,nds,nvs,NLINE,80) 


с Zone 6 Line 3 
NLINE=3 
DATA CIOS6,I,3),I21,8)/3,2,0,0,6,3,0,0/ 
write(80,97) (IC(6,I,NLINE),I=1,8) 


call zonelines(zone6x,zone6y,nds,nvs,NLINE,80) 
€ Zone 6 Line 4 

NLINE=4 

DATA (IC(6,I,4),I=1,8)/4,1,0,0,6,4,0,0/ 

write(80,97) (IC(6,I,NLINE),I=1,8) 


write(80,*) zoneb5x(nbot-1,nvs),zone5y(nbot-1,nvs) 
call zonelines(zone6x,zone6y,nds,nvs,NLINE,80) 
с 
С105Е (80) 
C kk koe oe oe o e oe oe oe oe oe o a ad oe oe oe oe ok ok oc oc ok oe oe oe oe oe oe oe oe ok oe oe oko oc oe oe oe ok ok ok ok ok ok ok oe oc oc oe oe ok ok oc oc oc oc oo a a a ad a ad ad ad a 
IN ORDER TO WRITE THE RESTART FILE, ZONES 2,3,5,6 NEED TO 
BE APPENDED WITH THE "N-1" COLUMN FROM THE UPSTREAM ARRAY 
SO THAT OVERLAPPING CONTINUITY CAN BE MAINTAINED 
Append zone 2 
CALL MAKEDUM(ZONE1X,ZONE1Y,NUS,NVS,ZONE2X,ZONE2Y,NTOP,NVS, 
-DZONE2X ,DZONE2Y , NTOPD,NVSD) 
C Append zone 3 
CALL MAKEDUM(ZONE2X,ZONE2Y,NTOP,NVS,ZONE3X, ZONE3Y,NDS,NVS, 
+DZONE3X , DZONE3Y ,NDSD ,NVSD) 
C Append zone 5 
CALL MAKEDUM(ZONE4X ,ZONE4Y,NUS,NVS,ZONESX ,ZONESY ,NBOT,NVS, 
+DZONESX , DZONESY , NBOTD, NVSD) 
C Append zone 6 
CALL MAKEDUM(ZONESX, ZONESY , NBOT, NVS, ZONE6X, ZONE6Y ,NDS,NVS, 
+DZONE6X , DZONE6Y ,NDSD,NVSD) 


оо а а 


C k k k k k k k k ak k k k k k k k k k k k k k k k k ak k k k k k ak k k ak k ak k ak ak ak ak ak k ak k k k ak ak k k ak ak k k ak k k k k k k k k k k k k 


c WRITE INMESH RESTART FILE 
IF (NUMIT.GT.1)THEN 
WRITE (*,’(A)’) "WRITING FILE: ^ // RESTFILE 
OPEN (UNIT=78,FILE=RESTFILE,FORM=*UNFORMATTED?, 
+ STATUS=’ UNKNOWN’ ) 
NZONE=6 
WRITE(78)NZONE 
E *x*xx*xxZONE 1 жжжжжжжж 
NZONE=1 
WRITE(78)NUS,NVS 
WRITE(78) ((ZONE1X(I,J),ZONE1Y(I,J),I-1,NUS),J21,NVS) 
WRITE(78) ((IC(NZONE,M,N),M=1,8),N=1,4) 


C 
С *x**xZONE 2*x*x*x*x***x 
NZONE-2 
WRITE(78)NTOPD, NVSD 
WRITE(78) ((DZONE2X(I,J),DZONE2Y(I,J),I-1,NTOPD) , J21,NVSD) 
WRITE(78) (CIC(NZONE,M,N) ,M=1,8) ,N=1,4) 
C 
E жжжж2( МЕ Зжжжжжжжж 
NZONE=3 
WRITE(78)NDSD,NVSD 
WRITE(78) ((DZONE3X(I,J),DZONE3Y(I,J),I=1,NDSD),J=1,NVSD) 
WRITE(78) ((IC(NZONE,M,N),M=1,8),N=1,4) 
С 


119 





****ZONE 4жжжжжжжж 
NZONE=4 
WRITE(78)NUS,NVS 
WRITE(78) ((ZONE4X(I,J),ZONE4Y(I,J),I=1,NUS) , Jz1,NVS) 
WRITE(78) (CIC(NZONE,M,N) ,M=1,8) ,N=1,4) 


жжжж20МЕ Бжжжжжжжж 
NZONE=5 
WRITE(78)NBOTD ,NVSD 
WRITE(78) ((DZONESX(I,J) ,DZONESY(I,J),I=1,NBOTD) , J=1,NVSD) 
WRITE(78) (CIC(NZONE,M,N) ,M=1,8) ,N=1,4) 


жжжж20МЕ бжжжжжжжж 
NZONE=6 
WRITE(78)NDSD,NVSD 
WRITE(78) ((DZONE6X(I,J),DZONE6Y(I,J),I-1,NDSD),Jz1,NVSD) 
WRITE(78) ((IC(NZONE,M,N),M=1,8),N=1,4) 


CLOSE(78) 
ENDIF 
ELSE 
write(*,*) ’INMESH input files have not been written.’ 
ENDIF 


(C 2k ak al ak ad Ak ak od ad a e al al ade al ak ae al al al ad e al al al a e ad ad ad a ak ak a ad ak de ae a ae a ae ae ae ae ae ae ae ae ae abe ae ae ae a k k kk k k k k 


WRITE(*,*) *FIT2D COMPLETE?” 
END 


С$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 
С$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 


о с СУС) 


сусу Су 


000000000 


This subroutine takes rans wake data and makes it in the defining 
line for the downstream grid zone boundary. 


SUBROUTINE RANSWAKE(DSL,XTE, YTE, 
XYIN,XOUT, YOUT , NOUT ,NUMIN) 

IMPLICIT NONE 

HEBPESDSLSXTE,YTE,PI,TEMP1,TEMP2 

INTEGER NOUT,NUMIN,NI,I,K 

PARAMETER (NI=200,PI=3.141592653589793E00 ) 


ARRAYS 


REAL XYIN(2,NI),XOUT(NI),YOUT(NI),YTEMP(NI) 
REAL XTEMP (NI) ,WAKCUB(4*NI) 


DEBUGGER WHAT IS COMING IN? 


WRITE(*,*) (XYIN(1,K),K=1,NUMIN) 

WRITE(*,*) ?xxx? 

WRITE(*, *) (XYIN(2,K) ,K=1,NUMIN) 

CHECK THE DATA: 
Does the stream line start aft of the trailing edge? 


IF (XYIN(1,1).LT.XTE)THEN 
игтъе(ж,ж) ? жжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжж? 


WRITE(*,*)'ERROR: RANS wake starts forward of TE’ 
ыгісе(ж,ж) жж жж ж acer Gio a aC iia? 
NUMIN=-1 

ENDIF 


Does the stream data end before the domain runs out? 


TEMP1=XTE+DSL 
IF (XYIN(1,NUMIN).GT.TEMP1) THEN 
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ооо ооо 


100 


ооо оо aa 


aa 


C 
C 


нгіте(ж, ж)  жжжжжжжжжжж жж жж жжжжжжжжжжжжжжжжжжжж? 
WRITE(*,*)'ERROR: RANS wake extends past flow domain’ 
Write(%* ,%*) II kK a ae? 
NUMIN--1 

ENDIF 


IF THE DATA IS OK THEN PROCEED THROUGH THIS: 
IF (NUMIN.GT.O)THEN 


FIRST THE LINE NEED TO BE FIXE UP A LITTLE. 

THE BEGINNING NEEDS TO START AT THE TRAILING EDGE 
OF THE FOIL AND THE END NEED TO BE AT THE DOWN 
STREAM LIMIT OF THE FOIL 


DO 100, I=2,NUMIN+1 
XTEMP(I)=XYIN(1,I-1) 
YTEMP (I)-XYIN(2,I-1) 

CONTINUE 


PUT BEGINNING AT THE TRAILING EDGE 
XTEMP(1)=XTE 
YTEMP(1)-YTE 


PUT THE END AT THE DOWNSTREAM LIMIT 
XTEMP (NUMIN+2)=TEMP1 
YTEMP (NUMIN+2)=YTEMP (NUMIN+1) 


SPLINE THE LINE AND GET THE CUBIC COEFFICIENTS 
CALL UGLYDK(NUMIN+2,1,1,XTEMP,YTEMP,0,0,WAKCUB) 


DEBUGGER 

WRITE(*,*) ’ONE OF THE CUBICS IS: ’,WAKCUB(5) 
WE WANT TO SEND NOUT X,Y POINTS BACK TO MAIN 
DSL IS HOW FAR TO MARCH TOTAL 
TEMP2=DSL/(FLOAT(NOUT)-1.0) 

DEBUGGER 

WRITE(*,*) 'INTERVAL IS: ’,TEMP2 

XOUT(1)=XTE 


DO 200 I-2,NOUT 

XOUT (I )=XOUT(I-1)+TEMP2 
DEBUGGER 

WRITE(*,*)'THE XOUT VALUE IS:’ ,XOUT(I) 
CONTINUE 


NOW GET THE CORRESPONDING Y VALUES 
CALL EVALDK(NUMIN-*2,NOUT,XTEMP,XOUT , YOUT, WAKCUB) 


WRITE(*,*)?xxx? 
WRITE(*,*)’***RANS DATA INCORPORATED IN TO GRID***? 
WRITE(*,*)?x*xx? 
WRITE(*,*)NUMIN,NOUT 
WRITE(*,*) ((XOUT(k) , YOUT(k)) ,K=1,NOUT) 
ENDIF 


RETURN 
END 


CF A te ыы ыы ыы 


C 
С 


This subroutine uses a two-d vortex lattice lifting line to 
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determine find the wake dividing line behind an arbitrary 
foil section. This subroutine is a generalization of the 
VLM2D code presented in the 13.04 course notes. 


SUBROUTINE ADAPT(PRESSX,PRESSY,SUCTX,SUCTY,WAKEX,WAKEY, 
+ NPRES , NSUCT,DSL , NWAKE , TUNPARAM) 

IMPLICIT NONE 

INTEGER NPRES,NSUCT,NWAKE,I,J,K,N,M,NI,NPTS,NPAN,IERR,NTOT 
INTEGER NIMAGEPR 

REAL PI 

PARAMETER (NI=200,PI=3.141592653589793E00) 

REAL DSL,DELS,LEN,DX,DY,TOP,sum,ENDX,U,V,TEMP,STEP 

REAL TUNPARAM, RTEMP, DXR, DYR, YIMAGE,signimage 


ARRAYS 


REAL PRESSX(NI),PRESSY(NI),SUCTX(NI),SUCTY(NI) 

REAL WAKEX(NI),WAKEY(NI),cambx(ni),camby(ni) 

REAL TOPCUBX(4*NI) ,BOTCUBX(4*NI) , CAMBCUBX (4*NI) 

REAL TOPCUBY (4*NI) , BOTCUBY(4*NI) , CAMBCUBY (4*NI) 
REAL WAKCUB(4*NI) 

REAL ARCTOP(NI),ARCBOT(NI),STEMP(NI),ARCAMB (NI) 

REAL TEMPTX(NI) , TEMPTY (NI) 

REAL TEMPBX (NI) , TEMPBY (NI) 

REAL SV(NI),SC(NI),DS(NI) 

REAL XV(NI),YV(NI) ,XC(NI) ,YC(NI) ,B(NI) , TX(NI) , TY(NI) 
REAL XN(NI),YN(NI),A(NI,NI) ,GAMMA(NI) , WKAREA(NI) ,dydx(ni) 
INTEGER IPIVOT(NI) 


What is the Y distance to the first image set? ) 
(the "Image Plane" is at half this distance 
YIMAGE-1/ (TUNPARAM) 


write(*,*)'the wall is at:' , ywall 
Spline Foil Offsets and determine cubic coefficients. 


Array ARC___ 1s returned with arclength params non-dimed from 0..1 


CALL PUGLYDK(NSUCT ,SUCTX , SUCTY , ARCTOP , TOPCUBX , TOPCUBY ) 
CALL PUGLYDK(NPRES , PRESSX , PRESSY , ARCBOT , BOTCUBX , BOTCUBY ) 


Extract Points along the top and bottom surface to find the 
mean camber line. 


HOW MANY POINTS TO EXTRACT? 
NPTS=40 


AT WHAT ARC COORDINATE LOCATIONS TO EXTRACT? 
write(*,*) arctop(nsuct) 
STEMP(1)=0.0 


DO 100 I=1,NPTS 
STEMP (1)=(FLOAT(I-1)/FLOAT(NPTS-1))*ARCTOP(NSUCT) 


write(*,*)stemp(i) 
CONTINUE 
EXTRACT THE VALUES ON THE SUCTION SIDE: 


CALL PEVALDK(NSUCT,NPTS ,STEMP , ARCTOP , TEMPTX , TEMPTY , TOPCUBX 
+, TOPCUBY ) 
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AT WHAT ААС COORDINATE LOCATIONS? 
STEMP(1)=0.0 


DO 200 I=1,NPTS 
STEMP(I)Z(FLOAT(I-1)/FLOAT(NPTS-1) ) *ARCBOT(NPRES) 

CONTINUE 

EXTRACT THE VALUES ON THE PRESSURE SIDE: 


CALL PEVALDK(NPRES ,NPTS ,STEMP , ARCBOT, TEMPBX , TEMPBY , BOTCUBX 


+, BOTCUBY) 


FIND THE MEAN CAMBER LINE(SIMPLE MEAN OF TWO COORDINATES): 


DO 300 I=1,NPTS 
CAMBX (I )=(TEMPTX (I )+TEMPBX (NPTS-I+1) )*0.50 
CAMBY (I )=(TEMPTY (I )+TEMPBY (NPTS-I+1) )*0.50 


debugger 
WRITE(88,*)CAMBX(I),CAMBY(I) 


300 CONTINUE 


debugger 


Spline the mean camber line to obtain the Cubic Coefficients 
Once again using the arc parameter scaling. 
CALL PUGLYDK(NPTS,CAMBX,CAMBY , ARCAMB, CAMBCUBX , CAMBCUBY) 


DETERMINE THE LOCATIONS OF THE VORTICES AND CONTROL POINTS 


How many vortices and control points should there be? 
A convergence study was performed by varying the number 
of panels. The difference in lift coefficient is less 
than 0.1/4 when stepping from 20 to 40 panels. So. 
20 panels is sufficient. 

NPAN-40 


How many pairs of image foils should there be? 
Obviously, the images should only be put in in pairs. 
This value will be adjusted until there is a converged 
lift coefficient. 


A convergence study was conducted. The relative error is 
approx 0.04% with 16 image pairs. This is sufficient 
for this application. 

NIMAGEPR=16 


Establish COSINE spacing on ARC Length Parameter 
this is straight out of the 13.04 notes. 

DELS is the COSINE spacing interval 
DELS=PI/NPAN 


DO 400 I=1,NPAN 
SV(I)20.50*(1.0-COS((I-0.50) *DELS)) 
SC(I)20.50*(1.0-COS(I*DELS)) 
DS(I)-PI*SQRT(SV(I)*(1.0-SV(I))) /NPAN 


400 CONTINUE 


EXTRACT THE X AND Y VALUES OF THE VORTICES AND CONTROL POINTS 


CALL PEVALDK(NPTS,NPAN,SV,ARCAMB, XV, YV, CAMBCUBX 


+,CAMBCUBY) 


CALL PEVALDK(NPTS,NPAN,SC,ARCAMB, XC, YC, CAMBCUBX 


+,CAMBCUBY) 
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To Panel the walls insert that stuff here: 
How far away are the walls:  TUNPARAM dictates 


How far do you want to panel up an down stream of 
the foil? DSL is already passed in. So a reasonalbe 
Lanse tom paneling mis LE - DSL to TE + DSL + 1 


What is a reasonable spacing: UNIFORM 


How many panels per WALL: 2*NPAN seems OK 


write locations out to datafile for inspection to fort.88 
debugger 
WRITE(88,*)'VARIABLES-XV,YV, XC, YC? 
WRITE(88,*) 'ZONE T-VORT' 
DO 500 I=1,NPAN 
WRITE(88,*)XV(I), YV(I),XC(I),YC(I) 
CONTINUE 


Calculate the UNIT NORMAL XN,YN at each control point 
unit normal not used for anything yet, but it may be 
useful in the future. 


Also, slope is required on the foil for V*n equation later 


write(88,*)'zone? 

DO 600 I=1,NPAN-1 
DX=XV(I+1)-XV(I) 
DYzYV(I-*1)-YV(I) 
LEN=SQRT (DX*DX+DY*DY ) 
XN(I)=-DY/LEN 
YN(I)=DX/LEN 
dydx(i)=dy/dx 

write(88,*) xn(i),yn(i) 

CONTINUE 


FIX UP THE VALUE AT THE TRAILING EDGE BECAUSE THERE IS NO 
VORTEX AFTER THE LAST CONTROL POINT. 


DX-XC(NPAN) -XV (NPAN) 
DY=YC(NPAN)-YV(NPAN) 
XN(NPAN)=-DY/SQRT (DX*DX+DY+*DY) 
YN(NPAN)=DX/SQRT(DX*DX+DY+*DY) 
dydx(npan)=dy/dx 
write(88,*)xn(npan),yn(npan) 


ALL THE UNIT NORMALS £ Slopes ARE NOW ESTABLISHED. 


COMPUTE THE INFLUENCE COEFFICIENTS A(N,M) AND THEN 
INVERT THE MATRIX. (linearized) 

The A matrix are the influence coefficients. Based on 
poor agreement with RANS, it was decided to go to the 
true locations of the vortices on the mean camber line 
vice a linearized surface. 


ТОР-1.0/(2.0жРІ) 

DO 700 М-1,МРАМ 

DO 700 M=1,NPAN 
rtemp=sgrt((xv(m)-xc(n))**2+(yv(m)-yc(n))**2) 
dyr=-(xc(n)-xv(m))/rtemp 
dxr=-(yc(n)-yv(m))/rtemp 
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A(N,M)-(TOP/(rtemp) ) *(dxr*xn(n)*dyr*yn(n)) 
write(89,*)a(n,m) 
The above include the influence of the vortex influence on 
the foil itself. Now, add in the images above and below 
the foils due to wall effects. Need Tunparam 


DO 710 K=1,NIMAGEPR 


Test K to see if 1t is a positive or negative image line 
if true then it is a positive image line 
if not true then it is a negative image 


if ((2*int(K/2)).eq.k) then 
signimage-1.0 
else 
signimage--1.0 
endif 
DEBUGGER 
if ((n.eq.1).and.(m.eq.1)) then 
write(*,*) signimage 
endif 
First Add in the Effect due to the image ABOVE (in pos y-dir) 


rtemp=sqrt((xv(m)-xc(n))*x*2+ 

% ((k*yimage* signimage*yv(m))-yc(n))**2) 
dyr--(xc(n)-xv(m))/rtemp 
dxr=-(yc(n)-(k*yimage+signimage*yv(m)))/rtemp 

A(N,M)= A(N,M)+signimage*(TOP/(rtemp))*(dxr*xn(n)+dyr*yn(n)) 
Next, Add in the Effect due to the image BELOW (in neg y-dir) 
rtemp-sqrt((xv(m)-xc(n))**2- 

y ((-k*yimaget signimage*yv(m))-yc(n))**2) 
dyr--(xc(n)-xv(m))/rtemp 
dxr=-(yc(n)-(-k*yimage+signimage*yv(m)))/rtemp 

A(N,M)= A(N,M)+signimage*(TOP/(rtemp))*(dxr*xn(n)+dyr*yn(n)) 


continue 


CONTINUE 
FACTOR THE A MATRIX 
CALL FACTOR(A,IPIVOT,WKAREA,NPAN,NI,IERR) 


COMPUTE THE w VELOCITY AT THE NTH CONTROL POINT 
THE B MATRIX OF A*GAMMA-B 


DO 800 N=1,NPAN 
B(N)=DYDX (N) 
write(*,*)dydx(n) 
CONTINUE 
NOW BACK SUBTITUTE TO EXTRACT THE VORTEX STRENGTHS 
CALL SUBST(A,B,GAMMA,IPIVOT,NPAN,NI) 


check the circulation distribution on the foil 


write out to fort.89 
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sum=0.0 
do 900 ı=1,npan 

sum-sum-tgamma(i) 
write(89,*)xv(i),gamma(i)/ds(i) 
continue 
write(*,*) ’Estimated CL=’ ,sum*2.0 
write(*,*) ’Using Vortex Lattice’ 


extract the wake 


The circulation is now known for each vortex point. 

To find the wake, "Step" off the trailing edge and 
evaluate the field point velocity. March off a distance 

in that direction. Evaluate the field point velocity again. 
adjust course. Keep a record of the path. That will be 
the line for the wake centerline. 


The following quantities are needed: 
NWAKE, DSL 
Want to calculate: WAKEX,WAKEY © NWAKE points between TE and DSL 


March down the wake first, spline the values and then extract 
the NWAKE points. 


Where is the Trailing edge? € pressx(1),pressy(y) 
Where do you want to stop? 0 pressx(1)+DSL £ corresponding y 
endx=pressx(1)+ds1+1.0 


WHAT STEP DO YOU WANT TO MARCH IN? (think of as scaled time) 
STEP=DSL/25.0 


START MARCHING 


I-1 
TX(1)=PRESSX(1) 
TY(1)=PRESSY(1) 

write(88,*) ’zone’ 

DO 910 WHILE (TX(I).LT.ENDX) 


FIND FIELD POINT VELOCITY DISTURBANCE DUE TO VORTICES 
REFERENCE NEWMAN PAGE 190 
U=1 to add in the freestream effect 
Ш-1 0 
у=0.0 
DO 920 N=1,NPAN 
TEMP-2*PI*((XV(N)-TX(I)) *(XV(N) -TXCI)) 4 YV(N) -TY(I))* 


^ (ТУСМ)-ТҮСІ?)) 


UZU-GAMMA (N) * (TY (I)-YV (N) )/ TEMP 
V=V-GAMMA(N) * (TX (I)-XV(N) ) /TEMP 


Add in the influence of the images 


DO 1950 K-1,NIMAGEPR 


Test K to see if it is a positive or negative image line 
if true then it is a positive image line 
if not true then it is a negative image 


if ((2*int(K/2)).eq.k) then 
signimage-1.0 

else 
signimage--1.0 
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endif 


a a 


DEBUGGER 
First Add in the Effect due to the image ABOVE (in pos y-dir) 


tempz2*PI*(((XV(N)-TX(I))**24 
((k*yimaget signimage*yv(n))-ty(i))**2)) 


U=U+signimage*GAMMA (N)*(TY(I)-(k*yimage+ 
signimage*yv(n)))/TEMP 
V-V-signimage*GAMMA(N)*(TX(I)-XV(N))/TEMP 


с Next, Add in the Effect due to the image BELOW (in neg y-dir) 
€ 
temp-2*PI*(((XV(N)-TX(I))**24 
^ ((-k*yimage-* signimagex*yv(n))-ty(i))**2)) 
U-U*signimage*GAMMA(N)* (TY(I)-(-k*yimage- 
2 signimage*yv(n)))/TEMP 
V=V-signimage*GAMMA (N) *(TX(I)-XV(N) )/TEMP 
1950 continue 
© 
с 
С 
€ 
€ 
920 continue 
C 
C CALCULATE THE NEXT LOCATION TO LOOK 
€ 
с Make the first step a baby step 
€ 


IF(I.LE.3)THEN 
STEP=STEP/8.0 
ELSE 
STEP=DSL/25.0 
ENDIF 


о 


Where is the next uake location? 


TX(I+1)=TX(I)+U*STEP 
TY(I+1)=TY(I)+V*STEP 
write(88,*)tx(i),ty(i),0,0 


I=I+1 
910 END DO 
NTOT=I 
e 
С SPLINE THE POINTS TO OBTAIN CUBIC COEFFICIENTS 
C 
C...SPLINE SPACING WITH FIXED SLOPE AT THE ENDS 
CALL UGLYDK(NTOT,1,1,TX,TY,0,0,WAKCUB) 
C...EVALUATE SPLINE TO FIND STEP SIZE AT INTERMEDIATE POINTS 
С 
G WHERE ARE THE OUTPUT POINTS? 
C 
WAKEX(1)=PRESSX(1) 
DO 950 I=1,NWAKE-1 
WAKEX(I+1)=WAKEX(I)+DSL/FLOAT((NWAKE-1)) 
950 CONTINUE 
С 
© EXTRACT THE POINTS FOR THE WAKEX AND WAKEY TO RETURN TO MAIN 
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CALL EVALDK(NTOT,NWAKE,TX,WAKEX , WAKEY , WAKCUB) 
write(88,*)’ZONE T=WAKEADS’ 

DO 960 I=1,10 

WRITE(88,*)WAKEX(I) ,WAKEY(I),0,0 

CONTINUE 


мгісе(ж,ж) ? жжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжюжжжжжжжжж) 
write(*,*) ’THE GRID BOUNDARIES HAVE BEEN ADAPTED TO THE WAKE’ 
write(*,*) ? USING VORTEX LATTICE METHOD’ 

Write(*,*) ? жжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжжж? 
RETURN 

END 


C HH PE) HE EE HR RA 
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SUBROUTINE FACTOR(W,IPIVOT,D,N,NDIM,IERR) 
IMPLICIT REAL(A-H,0-Z) 
DIMENSION W(NDIM,NDIM) , IPIVOT(*) ,D(*) 
IERR=0 
ВОТЬ М 
IPIVOT(I)=I 
ROWMAX=0. 
po М 
ROWMAX=MAX (ROWMAX , ABS(W(I,J))) 
CONTINUE 
IF (ROWMAX.EQ.0.) GO TO 999 
D(I)=ROWMAX 
CONTINUE 
NM1=N-1 
IF(NM1.EQ.0) RETURN 
DO 20 K=1,NM1 
J=K 
KP1=K+1 
IP=IPIVOT(K) 
COLMAX=ABS(W(IP,K))/D(IP) 
DO 11 I-KP1,N 
IP=IPIVOT(I) 
AWIKOV-ABS(W(IP,K))/D(IP) 
IF(AWIKOV.LE.COLMAX) GO TO 11 
COLMAX-AWIKOV 
Jg 
CONTINUE 
IF(COLMAX.EQ.O.) GO TO 999 
IPK=IPIVOT(J) 
IPIVOT(J)=IPIVOT(K) 
IPIVOT(K)=IPK 
DO 20 I-KP1,N 
IP=IPIVOT(I) 
W(IP,K)=W(IP,K)/W(IPK,K) 
RATIO=-W(IP,K) 
DO 20 J=KP1,N 
W(IP,J)-RATIO*W(IPK,J)-*W(IP, J) 
CONTINUE 
IF(W(IP,N).EQ.0.) GO TO 999 
RETURN 
IERR=2 
RETURN 
END 
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SUBROUTINE SUBST(W,B,X,IPIVOT,N,NDIM) 
IMPLICIT REAL(A-H,0-Z) 
DIMENSION W(NDIM,NDIM),B(*),X(*), IPIVOT(*) 
IF(NSGT.1) GO TO 10 
X(1)=B(1)/W(1,1) 
RETURN 
10 IP=IPIVOT(1) 
X(1)=B(IP) 
DO 15 K=2,N 
IP=IPIVOT(K) 
KMi-K-1 
SUM=0. 
DO 14 J=1,KM1 
SUM=W(IP,J)*X(J)+SUM 
14 CONTINUE 
X(K)=B(IP)-SUM 
15 CONTINUE 
X(N)=X(N)/W(IP,N) 
КМ 
DO 20 NP1MK=2,N 
KP 1=K 
K=K-1 
IP-IPIVOT(K) 
SUM=0.0 
DO 19 J=KP1,N 
SUM=W(IP,J)*X(J)+SUM 
19 CONTINUE 
X(K)=(X(K)-SUM)/W(IP,K) 
20 CONTINUE 
RETURN 


E A m. A O is e е 
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Krk 


*END OF SINGLE PRECISION VERSION OF DAVE GREELEY’S DIRECT SOLVER 
IR RRR RRR RRR RRR RRR RRR RRR RR OR RRR RRR RR FOR ЖжЖж Жж ЖЖЖЖЖЖЖ 


C 
HH HH 


с 
SUBROUTINE MAKEDUM(AX,AY,IMAXA, JMAXA, BX, BY, IMAXB, JMAXB, 


-XOUT , YOUT, IMAXOUT , JUAXOUT ) 


C 
C 
C This subroutine takes in two axially adjacent arrays and 
с appends the imaxa-1 vertical line of points to the left 
с hand side of the downstream array. 
с 
с Last Updated:11 November 
Е 
IMPLICIT NONE 
INTEGER IMAXA, JMAXA, IMAXB, JMAXB, IMAXOUT , JHAXOUT 
INWEGER I,NI,J 
PARAMETER(NI-200) 
C 
C ARRAYS 
C 
REAL AX(NI,NI),AY(NI,NI),BX(NI,NI),BY(NI,NI) 
REAL XOUT(NI,NI),YOUT(NI,NI) 
C 
C 
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300 
200 


DO 100 J=1,JMAXA 
XOUT(1,J)=AX(IMAXA-1,J) 
YOUT(1,J)=AY(IMAXA-1,J) 

CONTINUE 


DO 200 I=1,IMAXB 
DO 300 J=1,JMAXB 
XOUT(1+1,J)=BX(1,J) 
YOUT(I+1,J)=BY(I, J) 
CONTINUE 
CONTINUE 
IMAXOUT=IMAXB+1 
JMAXOUT=JMAXB 
RETURN 
END 


Са ааа ааа ааа аны ыы а-а 


00000000 


Q 00 


100 


300 


200 


400 


с 


SUBROUTINE ZONELINES(X,Y,IMAX,JMAX,NLINE,NFILE) 


This subroutine assists in writing output file for FIT2D 
to INMESH output file. 


X,Y are arrays containing points 
IMAX,JMAX are # of array elements 
NLINE is the line number in the zone 
NFIEL is the file output number 


IMPLICIT NONE 
INTEGER NI,NJ,I, IMAX,JMAX,J,NLINE,NFILE 
PARAMETER (NI=200,NJ=200) 

REAL X(NI,NJ),Y(NI,NJ) 

REAL X1(NI,NJ),X2(NI,NJ) , Y1(NI,NJ) , Y2(NI, NJ) 
REAL RI1,RI2,RJ1,RJ2,BETA 


IF(NLINE.EQ.1)THEN 

DO 100 J=1,JMAX 
WRITE(NFILE,*)X(1,J),Y(1,J) 

CONTINUE 

ENDIF 

IF(NLINE.EQ.3)THEN 

DO 300 J=1,JMAX 
WRITE(NFILE,*)X(IMAX,J) ,Y(IMAX,J) 

CONTINUE 

ENDIF 

IF(NLINE.EQ.2)THEN 

DO 200 I=1,IMAX 
WRITE(NFILE, *)X(I,1),Y(I,1) 

CONTINUE 

ENDIF 

IF(NLINE.EQ.4)THEN 

DO 400 I=1,IMAX 
WRITE(NFILE,*)X(I,JMAX) , Y(I, JMAX) 

CONTINUE 

ENDIF 


RETURN 
END 


С++ ЧАА АА АА 


C 


C 
C 


SUBROUTINE BORDERS(X,Y,X1,Y1,X2,Y2,X3,Y3,X4,Y4, IX,JY,NI) 


This subroutine takes in the x,y points which represent the four 
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осо со осо сос со ос о о суо соо о 


00000000 
a 


aaa 


осо со с с со со сс) 


a 


осо со со со о сс) 


lines surrounding a 2d box. It then writes them to a 2d array 
which represents the box boundaries. The "interior" of the 2d 
array will still be zeros. The interior is filled in another 
subroutine which interpolates between all the points. 


| Side 4 | 
| | 
| ZONE "n" | 
Side 1 | | Side 3 
| | 
| | 
| Side 2 | 
МЕ. с ee - | 
/N 
| Array reference to enter finite 
| interpolation scheme. 
JY | 
| > 


IMPLICIT NONE 
INTEGER I,IX,JY,NTEMP,NI 


ARRAYS 
REAL X(NI,NI),Y(NI,NI),X1(NI) , Y1C(NI) , X2(NI) , Y2(NI) 


REAL X3(NI) ,Y3(NI) ,X4(NI) , YA(NI) 
REAL X5(NI),Y5(NI),X6(NI),Y6(NI) 


Check line One: Determine if the order forward up or forward down. 


ntemp = O then count forward 
ntemp - 1 count backwards 


debugger 


write(*,*)'Debugging in borders y(jmax), y(1)' 
write(*,*) y1(jy),y1(1) 
NTEMP=O 
IF (Y1(JY).LT.Y1(1))THEN 
NTEMP=1 
ENDIF 


Now put line one in to x and y arrays 


DO 1609 Iz21,JY 

IF (NTEMP.EQ.O) THEN 
ЖШ І)-Х1С(1) 
Хе мат) 

ELSE 
X(1,I)=X1(JY-I+1) 
Y(1,1)=Y1(JY-1+1) 

ENDIF 


debugger 
if (i.eq.1)then 
write(80,*) "ZONE" 


endif 
write so.*) x(1,1),y(1,1) 
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a 


end debugger 


a 


100 CONTINUE 
Check line Two: Determine if the order is forward or backwards 


ntemp = O then count forward 

ntemp = 1 count backwards 

NTEMP=0 

IF (X2(IX).LT.X2(1))THEN 
NTEMP=1 

ENDIF 


оооо о 


а 


Now put line two in to x and y arrays 


DO 200 I=1, IX 

IF (NTEMP.EQ.0) THEN 
X(I,1)-X2(I) 
Y(I,1)=Y2(I) 

ELSE 
X(I,1)=X2(IX-I+1) 
Y(I,1)=Y2(IX-I+1) 

ENDIF 


debugger 
if (i.eq.1)then 
write(80,*) "ZONE" 
endif 
write(80,*) x(I,1),y(I,1) 


end debugger 


0000000000 


200 CONTINUE 
Check line Three:Determine if the order forward up or forward down. 


O then count forward 
1 count backwards 


ntemp 
ntemp 


0000 00 


NTEMP=0 

IF (Y3(JY).LT.Y3(1))THEN 
NTEMP=1 

ENDIF 


Now put line three in to x and y arrays 


о 


DO 300 I=1,JY 

IF (NTEMP.EQ.O) THEN 
X(ix,I)=X3(I) 
Ү(їх,1)=ҮЗ(1) 

ELSE 
X(ix,I)=X3(JY-I+1) 
Y(ix,I)=Y3(JY-I+1) 

ENDIF 


debugger 
if (i.eq.1)then 
write(80,*) "ZONE" 
endif 


write(80,*) x(ix,I),y(ix,I) 


end debugger 


0000000000 
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300 CONTINUE 
с 


схжжж 
с 
с 
€ Check line Four: Determine if the order is forward or backwards 
C 
C ntemp - O then count forward 
с ntemp = 1 count backwards 

NTEMP=O 

IF (X4(IX).LT.X4(1))THEN 

NTEMP=1 
ENDIF 


о 


Now put line four in to x and y arrays 


ЕШ 400 Izi.IX 

IF (NTEMP.EQ.0) THEN 
X(I,jy)=X4(I) 
Y(I,jy)=Y4(I) 

ELSE 
X(I,jy)=X4(IX-I+1) 
Y(I,jy)=Y4(IX-I+1) 

ENDIF 


debugger 
if (i.eq.1)then 
write(80,*) "ZONE" 
endif 


write(80,*) x(I,jy),y(I,jy) 


end debugger 


00000000400 


400 CONTINUE 


a 


CALL interp(X,Y,IX,JY) 
RETURN 
END 
Dr AA e Е тт в 
С 
SUBROUTINE interp(X,Y,IMAX,JMAX) 


C 
C 
с Below is a subroutine for doing ISOPARAMETRIC interpolation 
с for a single zone. It assumes that all four edges are 
€ already split up to their desired spacings. 
C This method is based on the simple isoparemetric method presented 
с in "FINITE ELEMENT PROCEDURES" by Bathe. 
с 
с 
с Written by John Dannecker for use in FIT2D.F 
с Last Modified: November 7, 1996 
© 
IMPLICIT NONE 
INTEGER NI,NJ,I, IMAX,JMAX,J 
PARAMETER (NI=200,NJ=200) 
REAL X(NI,NJ),Y(NI,NJ) 
REAL X1(NI,NJ),X2(NI,NJ),Y1(NI,NJ),Y2(NI,NJ) 
REAL RIL,RI2;RJ1,RJ2,beta 
c 
С 1-15 
с 125 
с DO ALL THE X(I,J) 
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200 


100 


aaa 


400 


300 


с 
с 


DO 100 I=2,IMAX-1 
RSE, 1)-X(1,1))/(XCIMAX, 1)-X(1,1)) 
ri2-(x(i,jmax)-x(1, jmax) )/ (x (imax, jmax)-x(1, jmax)) 

RI1ZRI1*(X(I,JMAX) -X(1, jmax) )/(XCIMAX, JMAX) -X(1, JMAX)) 


DO 200 J=2,JMAX-1 


betaz(Y(1,J)-Y(1,1))/(Y(1,JMAX)-Y(1,1)) 
RJ1-RJ1-(Y(IMAX,J)-Y(IMAX, 1) )/CY(IMAX,Jmax) -Y( IMAX, 1)) 
ее) 


CONTINUE 
DOSALL THE Y(1,J) 


DO 300 j=2,jMAX-1 
Rji=(y(1,j)-y(1,1))/(y(1,jmax)-y(1,1)) 
Rj2-(y(imax,j)-y(imax,1))/(y(IMAX, JMAX) -y Càmax, 1)) 


DO 400 i-2,iMAX-1 
beta-(x(i,i1)-x(1,1))/(x(imax,1)-x(1,1)) 
Ri2-Ri2*(x(i,jmax)-x(1, jmax))/ (x (IMAX, jmax)-x(1, jmax)) 
Y(I,J)=y(i,1)+(y(i,jmax)-y(i,1))*(rj1*(1.0-beta)+rj2*beta) 
CONTINUE 


CONTINUE 
RETURN 
END 


HH HH HH HH ++ 


с 


00000000 


Q0 000 ana ana 


aa 


SUBROUTINE FINDYPLUS (RESBL,RESWALL , REYNOLD , TUNPARAM , USL, DSL) 
This subroutine calculates the parameter y+ 

using the method in Anderson,Tannehill,Pletcher 
Computaional Fluid Mechanics and Heat Transfer, 1984. 

Y+ is checked against the user input values for cell height 
on the foil and wall. If the user cell spacing exceeds 

y* criteria, then cell spacing can be changed. 

IMPLICIT NONE 

REAL RESBL,RESWALL,REYNOLD, TUNPARAM ,USL,DSL , HEIGHT , NU 


REAL LF,LW, VEL, YPLUSF,REYWALL , YPLUSW, TEMP 
CHARACTER ANSWER*3 


MIT WATERTUNNEL CROSS SECTION HEIGHT (FEET) 
HEIGHT-20.0/12.0 

KINEMATIC VISCOSITY (NU FT^2/SEC) 

FRESH WATER Q0 70 DEGREES F.(SOURCE PNA 1967) 
NU=1.0519E-5 

CALCULATE CHORD LENGTH 

LF=TUNPARAM*HEIGHT 

Calculate wall length 


LW=TUNPARAM*HEIGHT#*(1.O+USL+DSL) 
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nA NAN 


C 


CALCULATE FLOW SPEED IN TUNNER (FT/SEC) 
VEL=NU*REYNOLD/LF 

CALCULATE YPLUS FOR FOIL 

YPLUSF=(3 .0/ (REYNOLD** (-0 .10) *SQRT (0.0227 )*VEL/NU) )/LF 
CALCULATE WALL BASED REYNOLD NUMBER 

REYWALL=VEL*LW/NU 

CALCULATE YPLUS AT WALL 
YPLUSW=(3.0/(REYWALL**(-0.10)*SQRT(0.0227)*VEL/NU))/LF 


COMPARE USER INPUT CELL HEIGHT TO Y+: 


debugger 


FORMAT(A12,F10.6,A6,F10.6) 
write(*,*) ’Cell Dim User Input Val Yplus (3)? 


write(*,99) 'On Foil: ',resbl,' ",yplusf 
write(*,99) ’On Wall: ’,reswall,’ ’,yplusw 
ON THE FOIL: 


TEMP=YPLUSF/RESBL 
IF (RESBL.GE.YPLUSF) THEN 
WRITE(*,*) ’USER INPUT CELL HEIGHT ON FOIL IS GREATER THAN’ 
WRITE(*,*) ’CALCULATED Y+(3) VALUE.’ 
WRITE(*,?'(A,$)?) 'DO YOU WANT TO CHANGE TO CALC Y+(3)?<n>:? 
READ(*,? (A)?) ANSWER 
IF ((ANSWER.EQ.?Y?).0R. (ANSWER. EQ. ^ y?^)) THEN 
RESBL=YPLUSF 
WRITE(*,*)'CELL HEIGHT ON FOIL SET TO Y+(3).’ 
ELSE 
WRITE(*,*)'CELL HEIGHT ON FOIL IS UNCHANGED. 
ENDIF 
ELSE 
WRITE(*,*)’?USER SPECIFIED CELL HEIGHT ON FOIL IS ADEQUATE’ 
ENDIF 


AT THE WALL 
TEMP=YPLUSW/RESWALL 
IF (RESWALL.GE.YPLUSW) THEN 
WRITE(*,*) USER INPUT CELL HEIGHT AT WALL IS GREATER THAN’ 
WRITE(*,*) ?’CALCULATED Y+(3) VALUE.’ 
WRITE(*,?’(A,$)’) ’DO YOU WANT TO CHANGE TO CALC Y+(3)?<n>:’ 
READ(*,’(A)’) ANSWER 
IF (CANSWER.EQ.’Y’).OR. (ANSWER.EQ.’y’)) THEN 
RESWALL=YPLUSW 
WRITE(*,*) ’CELL HEIGHT AT WALL SET TO Y+(3).? 
ELSE 
WRITE(*,*)'CELL HEIGHT AT WALL IS UNCHANGED." 
ENDIF 
ELSE 
WRITE(*,*)'USER SPECIFIED CELL HEIGHT AT WALL IS ADEQUATE’ 


ENDIF 


RETURN 
END 


С+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 
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SUBROUTINE TUND(XBEG, YBEG , XEND , YEND , RESBL , RESWALL , XOUT , YOUT , 
*NPOINTS,PACKING) 


LAST MODIFIED:29 OCTOBER 1996 


This subroutine generates the vertical spacing of points on 
lines above and below the foil. 


000000 


IMPLICIT NONE 

REAL PI,TWOPI,ZERO,ONE,HALF 

INTEGER NI,I,NNEARBL,NNEARWALL , NRREMAIN 

PARAMETER (PI=3.14159,TWOPI=6.28318, NI=200, ZERO=0.0, ONE=1.0) 
INTEGER NPOINTS 

REAL XBEG, XEND, YBEG, YEND,RESBL,RESWALL, PACKING 

REAL DELX,DELY,TEMPX,TEMPY,DLEFT,DRIGHT 

real xvect,yvect,hypo 


€) CO 


ARRAYS 
REAL XOUT(NI),YOUT(NI),XLIN(NI),YLIN(NI),XREM(NI) , YREM(NI) 


Determine spacing of points for boundary layer along foil 


How many points close packed near foil and at the wall? 


о 0 O C) 


NNEARBL=INT(PACKING*FLOAT(NPOINTS)) 
NNEARWALL=NNEARBL 
NREMA IN=NPOINTS-NNEARWALL-NNEARBL+2 


Size the cells on the foil: 


aa 


XOUT(1)=XBEG 

YOUT(1)-YBEG 

DELX=XEND-XBEG 

DELY-YEND-YBEG 

hypo=sqrt (delx*delx+dely*dely) 
xvect=delx/hypo 

yvect=dely/hypo 

debugger 
write(*,*)xvect,yvect,delx,xend,xbeg 


Each cell is 25% larger than the previous cell. 


000000 


DO 100 I=1,NNEARBL-1 
XOUT(I-*1)2XOUT(I)*(1.25**(I-1))*RESBL*xvect 
YOUT(I+1)=YOUT(I)+(1.25**(I-1))*+RESBL*yvect 

с write(*,*)i+1,XOUT(I+1),YOUT(I+1) 
100 CONTINUE 
DLEFT=(1.25**(I-1))*RESBL 


E Size the cells at the wall: 
| XOUT (NPOINTS)-XEND 
YOUT(NPOINTS)=YEND 
E NOTE: This is a REVERSE counter!!!!! 
: DO 200 I-1,NNEARWALL- 1 
E Each cell is 25% larger than the previous cell. 


XOUT (NPOINTS-I )=XOUT (NPOINTS-I+1)-(1.25**(I-1) ) *RESWALL*xvect 
YOUT(NPOINTS-I)-YOUT(NPOINTS-I*1)-(1.25**(I-1))*RESWALL*yvect 
c write(*,*)npoints-i,XOUT(NPOINTS-I), YOUT(NPOINTS-I) 
200 CONTINUE 
DRIGHT2(1.25**(I-1)) *RESWALL 
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Send the remainder of the line to FNSPLT: 


First need to generate a small line 


Endpoints of line: XOUT(NNEARBL) , YOUT(NNEARBL) 
XOUT (NPOINTS-NNEARWALL) , YOUT(NPOINTS-NNEARWALL) 
FNSPLT WILL NEED 5 POINTS 
TEMPX-XOUT (NPOINTS-NNEARWALL--1)-XOUT (NNEARBL) 
TEMPY=YOUT (NPOINTS-NNEARWALL+1)-YOUT(NNEARBL) 


debugger 
write(*,*)tempx,tempy , XOUT(NPOINTS-NNEARWALL) 


XLIN(1)=XOUT(NNEARBL) 

YLIN(1)=YOUT(NNEARBL) 

DO 300 I=1,4 
XLIN(I+1)=XLIN(I)+0.25*TEMPX 
YLIN(I+1)=YLIN(I)+0.25*TEMPY 

write(*,*) ’Straight Line’ 
write(*,*)xlin(i),ylin(i) 

CONTINUE 


CALL FNSPLT(5,NREMAIN,DLEFT,DRIGHT,XLIN,YLIN,XREM,YREM) 


FILL OUT THE ARRAY 

DO 400 I-1,NREMAIN 
XOUT(NNEARBL+I-1)=XREM(I) 
YOUT(NNEARBL+I-1)=YREM(I) 

CONTINUE 


debugger 
WRITE(*,*)NNEARBL,NNEARWALL , NPOINTS,NREMA IN 


debugger 

write(79,*) ?20МЕ? 

do 500 i-1,npoints 
write(79,*) xout(i),yout(i) 
continue 

end debugger 


RETURN 
END 


Coon bee ebbe eee ee ie eed e een ed e e eR b ee e rn 


(0020200000 00400 


SUBROUTINE TUNC(XBEG, XEND , YVALL, NPOINTS , XIN , YIN, XOUT , YOUT) 
LAST MODIFIED:08 NOVEMBER 1996 
This subroutine takes the cell spacing along the upstream line 


in the center of the tunnel, the foil points and the down 
stream line and spaces them along tunnel walls. 


XBEG = BEGINNING XCOORD OF LINE 
XEND = ENDING XCOORD OF LINE 


YWALL = Y VALUE AT WALL 
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NPOINTS - NUMBER OF ARRAY VALUES COMING IN 
XIN,YIN, - VALUES COMING IN THAT WILL BE PROJECTED 


XOUT,YOUT- X AND Y VALUES FOR THE LINE ON THE WALL 


IMPLICIT NONE 

REAL PI,TWOPI,ZERO,ONE,HALF 

INTEGER NI,I 

PARAMETER (PI=3.14159,TWOPI=6.28318, NI=200, ZERO=0.0, ONE=1.0) 
INTEGER NPOINTS 

REAL XEND, XBEG,DELX,DELY,TEMP3,SLENGTH , TEMP4, YWALL 

REAL XIN(NI), YIN(NI), XOUT(NI), YOUT(NI),DELARC(NI) 


Compute arc length of the input array and arc length between 
points. 


SLENGTH-ZERO 

DELX=ZERO 

DELY=ZERO 

TEMP3=ZERO 

DO 100 I-1,NPOINTS-1 
DELX=XIN(I+1)-XIN(I) 
DELY=YIN(I+1)-YIN(I) 
TEMPS-DELX*DELX-*DELY*DELY 
DELARC(I)=SQRT(TEMP3) 
SLENGTH=SLENGTH+DELARC(I) 

CONTINUE 


debugger 
temp4=xin(npoints)-xin(1) 
write(*,*)slength,temp4,ywall,XBEG,XEND 


do algebraic translation of point spacing 


TEMP4=ZERO 

XOUT (1)=XBEG 

YOUT(1)=YWALL 

DO 200 I=1,NPOINTS-1 
TEMP4=DELARC(I)/SLENGTH 
XOUT(I+1)=X0OUT(I)+TEMP4*(XEND-XBEG) 
YOUT(I+1)=YWALL 

CONTINUE 


debugger 


write(79,*) ’ZONE’ 

do 300 i=1,npoints 
write(79,*) xout(i),yout(i) 
continue 


end debugger 


RETURN 
END 


Сани 


000000 


SUBROUTINE WALLFIND(SCALE,YUP,YBOT) 
LAST MODIFIED: 29 OCTOBER 96 


This subroutine finds the upper and lower y coordinates of the 
tunnel walls 
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IMPLICIT NONE 

REAL PI,TWOPI,ZERO,ONE,HALF,TWO 

INTEGER NI 

PARAMETER (PI=3.14159,TWOPI=6.28318,NI=10, 
+ZERO=0.0,ONE=1.0,TWO=2.0) 

REAL SCALE, YUP, YBOT 


YUP=ONE/ (SCALE*TWO) 
YBOT=-ONE/ (SCALE*TWO) 


debugger 


write(*,*) yup,ybot 


000000 


RETURN 
END 
Ire b e p bee т АА Аа Б Е ОР А 
SUBROUTINE TUNB(XL,XLOC,YLOC,PLINEX,PLINEY,NP) 
LAST MODIFIED: 24 OCTOBER 96 


This subroutine defines the upstream or downstream lines extending 
from the foil leading or trailing edge. 


XL = LENGTH OF LINE 
XLOC,YLOC = BEGINNING OF LINE 


PLINEX,PLINEY - OUTPUT X AND Y 


аза асса ао саса 


NP - NUMBER OF POINTS OUT 
IMPLICIT NONE 

REAL PI,TWOPI,ZERO,ONE,HALF 

INTEGER NP 

PARAMETER (PI=3.14159,TWOPI=6.28318, ZERO=0.0, ONE=1.0) 
REAL XL,XLOC,YLOC,TEMP 

INTEGER I 

REAL PLINEX(NP),PLINEY(NP) 


debugger 
write(*,*) Майе іт іп to TUNB, LIMIT IS:', XL 


о о о сус) 


MAKE IT HORIZONTAL 
DO 100 I=1,NP 
PLINEY(I)-YLOC 
100 CONTINUE 


debugger 
write(*,*) xloc,yloc 


ооо о со 


The endpoints of the line are xloc,yloc and xloc+temp,yloc 
IF (XLOC.LE. ZERO) THEN 

PLINEX (NP )=XLOC-XL 

PLINEX(1)=XLOC 
ENDIF 


IF (XLOC.GT.ZERO)THEN 
PLINEX (NP )=XLOC+XL 
PLINEX(1)-XLOC 

ENDIF 


aa 


Generate some points between the two endpoints. 
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000000 


TEMP=(PLINEX(NP)-PLINEX(1))/FLOAT(NP-1) 

DO 200 I=1,NP-2 
PLINEX(I+1)=PLINEX(I)+TEMP 

CONTINUE 


DEBUGGER 
write(79,*) 'ZONE' 
do 5020 i=1,NP 
write(79,*) PLINEX(i),PLINEY(i) 


5020 continue 


RETURN 
END 


HH HH ++ 
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SUBROUTINE ARCJOINER(N1,N2,NOUT,X1,X2,Y1,Y2,XOUT , YOUT) 
last modified: 25 October 96 


This subroutine takes in two arcs and joins them end to end 

in to one array. 

IMPLICIT NONE 

REAL PI,TWOPI,ZERO,ONE,HALF 

INTEGER NI,I 

PARAMETER (PI=3.14159, TWOPI=6.28318, NI=200, ZERO=0.0, ONE=1.0) 
INTEGER Ni, N2,NOUT 

REAL X1(NI), Yi1(NI), XOUT(NI), YOUT(NI), X2(NI),Y2(NI) 


debugger 
write(*,*) ’Made it in to ARCJOINER:’, ni,n2 


DOLIO I=1,N1 
XOUT(I)=X1(I) 
YOUT(I)=Y1(I) 

CONTINUE 

DO 20 I=2,N2 
XOUT(I+N1-1)=X2(I) 
YOUT(I+N1-1)=Y2(I) 

CONTINUE 


funny counting here elmininates overlapping data points 
NOUT=N1+N2-1 
debugger 


write(*,*) "Points Out',nout 


RETURN 
END 


Са а-а ааа ана ыы а-а ана а-а а-а-а ыы ыы ааа 


Сз СУ СУ СУ СУ СУ СУ СУ 


SUBROUTINE ARCMAKER(X,Y,NUM1,NUM2, XO, YO) 

last modified: 29Sept96 
THIS SUBROUTINE TAKES IN LONG ARRAYS X, AND Y, 
AND RETURNS PORTIONS OF X AND Y SPECIFIED BY 
INTEGERS NUM1 AND NUM2. ХО AND YO ARE FILLED 
BY THE FIRST NUM1-NUM2 ELEMENTS 


IMPLICIT NONE 
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REAL PI,TWOPI,ZERO,ONE,HALF 

INTEGER NI,I 

PARAMETER (PI=3.14159,TWOPI=6.28318, NI=200, ZERO=0.0, ONE=1.0) 
INTEGER NUM1,NUM2 

REAL X(NI), Y(NI),XO(NI), YO(NI) 


debugger 
write(*,*)NUM1,NUM2 

DO 10 I=1,(NUM2-NUM1+1) 
XO(I)=X(NUM1+I-1) 
YO(I)=Y(NUMi+I-1) 


debugger 
write(*,*) XO(I),YO(I) 


CONTINUE 
RETURN 
END 


С+++++++++++++++++++++++++++++++++++++++-+++++++++++++++++++++++++++++++ 
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SUBROUTINE FOILZONE(XI,NPOINTS,MIDPRES,MIDSUCT,NOSE) 
last modified: 29 Sept 96 


THIS SUBROUTINE LOOKS AT ALL THE FOIL GEOMETRY POINTS. IT SEARCHES 
FOR THE LEADING EDGE OF THE FOIL AT THE SET AOA. IT SEARCHES FOR 
POINTS ON THE SUCTION AND PRESSURE SIDES THAT ARE CLOSEST TO THE 
MIDCHORD POINT. ONCE THESE POINTS ARE IDENTIFIED, THE INTEGER 

OF THE ARRAY POSITION OF THESE POINTS IS RETURNED TO MAIN PROG 


MIDSUCT 
NOSE < ЕЕ ЕЕ |) (foil schematic) 


MIDPRES 


IMPLICIT NONE 

REAL PI,TWOPI,ZERO,ONE,HALF 

INTEGER NI,I 

PARAMETER (PI=3.14159,TWOPI=6.28318, NI=200, ZERO=0.0, ONE=1.0) 
INTEGER NPOINTS,MIDPRES , MIDSUCT, NOSE, TEMP 

REAL CHORDMID 

REAL XI(NI) 


XI,YI ARE FOIL OFFSETS 
NPOINTS : TOTAL NUMBER OF FOIL OFFSETS 
MIDPRES ,MIDSUCT AND NOSE AS INDICATED ON SCHEMATIC 


FIND THE LEADING EDGE: 


TEMP=0 
DO 10 I=1,(NPOINTS-1) 
IF ((XI(I)).GE.(XI(I+1))) THEN 
TEMP=I+1 
ENDIF 
CONTINUE 
NOSE=TEMP 


FIND THE MIDDLE OF THE PRESSURE SIDE 


TEMP=0 
CHORDMID=((XI(NOSE)+XI(1)))/(2.0) 


141 





27% 


осо со со 


20 


AN 


30 


00000 


C 


debugger 
write(*,*) ’pressmid=’, chordmid 


DO 20 I=1,(NOSE-1) 
IF (XI(I).GT.CHORDMID) THEN 
TEMP-I 
ENDIF 
CONTINUE 
MIDPRES-TEMP 


FIND THE MIDDLE OF THE SUCTION SIDE 


TEMP=0 
CHORDMID=(XI(NOSE)+XI(NPOINTS))/(2.0) 
write(*,*) ’suctmid=’, chordmid 
DO 30 I=NOSE, (NPOINTS-1) 

IF (XI(I).LT.CHORDMID) THEN 

TEMP-I : 

ENDIF 
CONTINUE 
MIDSUCT=TEMP 


debugger 
WRITE(*,*) NPOINTS ,NOSE,MIDSUCT,MIDPRES 


RETURN 
END 


Ся 4 ЧЕН ЧЕЧЕНА НЧЕ + + + + + 


Q00000000055000»050»0 


ann 


C 


SUBROUTINE NACACONV(N,X,Y,T,R) 
last modified: 23 Sept 96 
THIS SUBROUTINE CONVERTS NACA FOIL GEOMETRY TO X,Y GEOMETRY 


NUMBER OF STATIONS 
STATION 

CAMBER 

THICKNESS 

L.E. RADIUS 


N 
X 
Y 
T 
R 
RETURNS 
N = NUMBER OF POINTS (.GT. 2*STATIONS) 
X,Y = GEOMETRY COORDINATES 
IMPLICIT NONE 
REAL PI,TWOPI,ZERO,ONE,HALF 
INTEGER NI 
PARAMETER(PI=3.14159,TWOPI=6.28318, NI=200) 
INTEGER I,N 
REAL R 
REAL X(NI),Y(NI),T(NI) 


THIS SUBROUTINE NOT YET FULLY IMPLEMENTED 
WRITE(*,*) ' THE NACA CONVERSION ROUTINE IS NOT YET IMPLEMENTED’ 


WRITE(*,*) ’ Pretty lame, huh?’ 
END 


СНЕ HHHH 


С 
C 


SUBROUTINE ROTANGLE(N,X,Y , AOA, PIVOTX,PIVOTY) 


last modified: 26 Sept 96 
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anaaaaaaaaaaaN 


aaa 


101 


C 


THIS SUBROUTINE TAKES X,Y FOIL DATAPOINTS AND ROTATES THEM ABOUT 
A PIVOT POINT TO SET FOIL TO PRESCRIBED ANGLE OF ATTACK(AOA) 


N = NUMBER OF POINTS DEFINING FOIL 
X = X COORDINATE OF POINT (X/C) 
Y = Y COORDINATE OF POINT (Y/C) 


AOA = FOIN ANGLE OF ATTACK IN DEGREES REF TO NOSE TAIL LINE 
PIVOT = X/C OF ROTATION POINT OF FOIL 


XL,YL ARE REF TO PIVOT POINT 
R IS RADIUS TO PIVOT POINT 


IMPLICIT NONE 
REAL PI,TWOPI,ZERO,ONE,HALF 
INTEGER NI 
PARAMETER(PI=3.14159,TWOPI=6.28318, NI=200) 
INTEGER I,N 
REAL AOA, PIVOTX, PIVOTY, XL, YL, R, BETA 
REAL X(NI),Y(NI) 


WRITE(*,*) " FOIL BEING ROTATED TO ANGLE OF ATTACK:' , AOA 
WRITE(*,*) ? PIVOT POINT IS X COORD © (X/C):’ , PIVOTX 
WRITE(*,*) ’ PIVOT POINT IS Y COORD © (Y/C):’ , PIVOTY 


ROTATE FOIL POINT BY POINT 


DO 101 I =1,N 
XL=X(I)-PIVOTX 
YL=Y(I)-PIVOTY 
R = SQRT(XL*XL+YL*YL) 
BETA = ATAN2D(YL, XL) 


X(I) = R*COSD(BETA-AOA) 
Y(I) = R*SIND(BETA-AOA) 
CONTINUE 

RETURN 

END 


Catt ttt ett HH HH HH HH + 


C 


c 
c 


10 


c 
c 


SUBROUTINE NORMFOIL(L,N,X,Y) 


LAST MODIFIED: 24 OCTOBER 96 
IMPLICIT NONE 
REAL PI,TWOPI,ZERO,ONE,HALF 
INTEGER NI,I 
PARAMETER (PI=3.14159,TWOPI=6.28318, NI=200, ZERO=0.0, ONE=1.0) 
INTEGER N 
REAL L 
REAL X(NI),Y(NI) 
DO 107 I=1,N 
X(I)-X(I)/L 
Y(I)=Y(I)/L 
CONTINUE 
RETURN 
END 


HH HH HH + 


C 
C 
C 


SUBROUTINE FNSPLT(NIN,NOUT,DS1,DS2,XI,YI,XO, YO) 


ӨЫ 2747/84 5% 5 9 8 я © 9 е в €* в ә е өӨ 6 е ва € € 8 € € * « ес ео ө эе а ә €9 єс 8 © € 6 ө оо о о ае обоо е «e s 9 e оо 


THIS PROGRAM CONVERTED IN TO A SUBROUTINE IN "FIT2D.F" ! 


143 





! BY JOHN DANNECKER. THIS PROGRAM IS ORIGINALLY WRITTEN е 
BY S.D. BLACK, MIT, WITH OTHER CREDITS AS INDICATED ! 


е 
> о хо особо 9 16 ө е 4 4 е 96 4 ө 5 ө v ө а: ” 4 68 9:9 е ғ 9 ә е e 6 “4 © © «6 « 6 » 5 э % 


ОО ОО 


FANCY SPLITTING PROGRAM FOR DIVIDING UP A CURVE 
THERE ARE TWO METHODS THIS PROGRAM CAN BE RUN 


INPUT 1: 
DATA SET (X,Y) FOR THE CURVE 
DS1 - INTERVAL AT LEFT HAND END POINT 
DS2 - RIGHT HAND END POINT INTERVAL 
SET NOUT « O 


INPUT 2: 

SAME AS INPUT 1 BUT WITH THE NUMBER OF POINTS SPECIFIED 
OUTPUT 

X,Y OF SPLIT UP CURVE BETWEEN THESE INTERVALS 


PURPOSE 
USEFUL FOR CUSTOM GRIDDING OF CURVES 


THE SLOPE OF THE INTERVAL DS# IS HELD TO BE ZERO AT THE END POINTS 
THE PROGRAM ASSUMES THE ENDS OF THE CURVE ARE AT THE ENDS OF THE INPUT 


WRITTEN 1/8/95 S. BLACK, MOD BY J. DANNECKER 9/28/96 


This is the NEW and Improved FNSPLT modified by 
S. Black 25 October 1996 


oo oe onen nn na nanaaanaaanan 


PARAMETER (NI=200) 

REAL XI(NI),YI(NI),XO(NI),YO(NI),CUB(4*(NI-1)) 
REAL SI(NI),CUB1(4*(NI-1)),CUB2(4*(NI-1)) 
REAL A(3),B(3),C(NI),D(NI),E(NI) 

CHARACTER FIN*20 

CHARACTER PROMPT2*34 


C 
C 
с 
с debugger 
e write(*,*)'made it in to fnsplt', nin,nout 
СИНИЕ НИЕ НИИ ЕН HE SES 
с 
C CALCULATE LENGTH OF CURVE BASED ON INPUT POINTS 
SL=0.C 
DO 10 I-2,NIN 
SL-SL*SQRT((XI(I)-XI(I-1))**2* (YI(I)-YI(I-1))**2) 
10 CONTINUE 
C CALCULATE THE NUMBER OF POINTS REQ’D IF NOT SPECIFIED 
IF (NOUT.LE.O) THEN 
DSAVG=(DS1+DS2)/2.0 
NOUT=INT(SL/DSAVG)+1 
END IF 
SPLINE INPUT ARRAY PARAMETRICALLY (BOTH X AND Y) 
Array SI is returned with arclength parameters non-dimed from 0..1 
CALL PUGLYDK(NIN,XI,YI,SI,CUB1,CUB2) 
SET UP ARRAYS CONTAINING SPACING 
Array A Contains pointers for the first, middle and last 
steps. 
Array B M rins the dS values at the two ends and a guess 
at what the step size in the middle should be. 
Array C contains integers cooresponding to the NOUT-1 steps 


aa 


Са (а со ка Ка ср с) 
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pe 
| ы "s. 6 
nh 
1 

1d 


ГЕ 


n 


[^ 


C 


C 
C 
C 
C 
C 
C 
C 


C 
C 


C 
C 


C 


C 
с 
С 


DSAVG=(SL-DS1-DS2)/FLOAT(NOUT-3) 
A(1)=FLOAT(1) 
A(2)-FLOAT(NOUT)/2.0 
A(3)=FLOAT(NOUT-1) 
В(1)=р51 
B(2)=DSAVG 
B(3)=DS2 
DO 20 J=1,NOUT-1 
C(J)=FLOAT(J) 
20 CONTINUE 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
DSMIN=0.1*MIN(DS1,DS2) 
...CONVERGENCE LOOP FOR APPROPRIATE VALUE OF B(2) 

The step size is splined with the slope held constant at the 
ends. There must be NOUT-1 steps, with values of DS1 and 
DS2 at the two ends. The mid-range step size is varied until 
the sum of all the steps equals the arc length desired (SL). 
By fixing the slope at the ends the program tries to ensure 
that step sizes don't increase too rapidly. 

DO 100 KK=1,1000 
IERR=O 
...SPLINE SPACING WITH FIXED SLOPE AT THE ENDS 
CALL UGLYDK(3,0,0,A,B,0,0,CUB) 
. . EVALUATE SPLINE TO FIND STEP SIZE AT INTERMEDIATE POINTS 
CALL EVALDK(3,NOUT-1,A,C,D,CUB) 
. . . CALCULATE LENGTH COVERED BY NEW SET OF STEPS 
SLC=0.0 
DO 40 I-1,NOUT-1 
SLC=SLC+D(I) 
IF (D(I).LT.O.O) IERR=1 
40 CONTINUE 
...IF LENGTH IS NOT CLOSE TO TOTAL LENGTH NEEDED, SHIFT B(2) 
....AND REITERATE 
IF (IERR.EQ.1) THEN 
B(2)=0.0 
ELSE IF (ABS(SLC-SL).GT.DSMIN) THEN 
DIFF=SLC-SL 
B(2)=B(2)-DIFF/NOUT 
ELSE 
GOTO 101 
END IF 
100 CONTINUE 
WRITE(*,*)?***WARNING*** FNSPLT DID NOT CONVERGE’ 
WRITE(*,*)’ ABS(SLC - SL) > TOL’ 
WRITE(*,*)’ SLC, SL, TOL = ’,SLC,SL,DSMIN 
WRITE(*,*)NIN,NOUT,DS1,DS2 
IF (SLC.GT.SL) THEN 
WRITE(*,*)'Recommend using fewer points or smaller dS values’ 
ELSE 
WRITE(*,*)'Recommend using more points of larger dS values’ 
END IF 
101 CONTINUE 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
...SHIFT LENGTH TO EXTEND FROM 0.0 TO 1.0 
E(1)20.0 
DO 50 I-2,NOUT-1 
E(I)=D(I-1)/SLC+E(I-1) 
50 CONTINUE 
E(NOUT)=1.0 
...EVALUATE OUTPUT POINTS AT CALCULATED INTERVAL 
CALL PEVALDK(NIN,NOUT,E,SI,XO,YO,CUB1,CUB2) 


ааатшаттшынынытышышнытынынаыныты 
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RETURN 


END 
С------------------------------------------------------- С 
SUBROUTINE DRIVDK(NIN,NOUT,XIN,XOUT,DYDX,D2YDX,A) 
C APRIL 1975 SPLINE PROGRAM SERIES J.E.KERWIN 
REAL XIN(*),XOUT(*) ,DYDX(*) ,D2YDX(*) , A(*) 
NMi=NIN-1 
J=1 


DO 3 N=1,NOUT 
IF(XOUT(N).GE.XIN(2)) GO TO 4 
Тай 
GO TO 5 
4 IF(XOUT(N).LT.XIN(NM1)) GO TO 6 
J=NM1 
GO TO 5 
6 IF(XOUT(N).GE.XIN(J+1)) GO TO 7 
5 H1=X0UT(N)-XIN(J) 
H2-H1**2 
J2=J+NM1 
J3=J2+NM1 
DYDX(N)23.0*A(J) *H242.0*A(J2) *H1-A(J3) 
D2YDX(N)26.0*A(J) *H142.0*A(J2) 
GO TO 3 
ү J=J+1 
GO TO 6 
3 CONTINUE 
RETURN 
END 
EENCCCESEEEBEECCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCOCCC 
SUBROUTINE EVALDK(NIN,NOUT,XIN,XOUT,YOUT,A) 
C APRIL 1975 SPLINE PROGRAM SERIES J.E.KERWIN 
REAL XIN(*),XOUT(*) , YOUT (*) , A(*) 
NM1=NIN-1 
MOUT=IABS (NOUT) 
IF(NOUT.GT.O) GO TO 1 
DEL-(XIN(NIN)-XIN(1))/(MOUT-1) 
DO 2 N=1,MOUT 
2 XOUT(N)=XIN(1)+(N-1)*DEL 
1 Jai 
DO 3 N=1,MOUT 
IF(XOUT(N).GE.XIN(2)) GO TO 4 
и 
GO TO 5 
4 IF(XOUT(N).LT.XIN(NM1)) GO TO 6 
J=NM1 
GO TO 5 
6 IF(XOUT(N).GE.XIN(J+1)) GO TO 7 
9 IF (XOUT(N).LT.XIN(J)) GO TO 8 
5 H1=XOUT(N)-XIN(J) 
H2=H1**2 
H3=H1*H2 
J2=J+NM1 
J3=J2+NM1 
J4=J3+NM1 
YOUT(N)=A(J)*H3+A(J2)*H2+A(J3)*H1+A(J4) 
GO TO 3 
7 J=J+1 
GO TO 6 
8 J=J-1 
GO TO 9 
3 CONTINUE 
RETURN 
END 
C 


С+++++++++++++++++++++++++++++++++++++++++++++++++++++++-+++++++++++++++++ 
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C 


SUBROUTINE INTEDK(NIN,XIN,XL,XU, YDX, XYDX, XXYDX , A) 
AUGUST 1975 SPLINE PROGRAM SERIES S.-K.TSAO 
REAL XIN(*),A(*) 

NM1=NIN-1 

BMSCXESPT.XIN(1)) GO TO 2 
DO 1 N=1,NIN 
IF(XL.GE.XIN(N)) GO TO 1 
JL=N-1 

GO TO 3 

CONTINUE 

JL=NM1 

GO TO 3 

Sd 

IF(XU.GE.XIN(NIN)) GO TO 4 
Jet 

IF(XU.LE.XIN(1)) GO TO 6 
DO 5 N-JL,NIN 
IF(XU.GE.XIN(N)) GO TO 5 
JU=N-1 

GO TO 6 

CONTINUE 


Hi=XL-XINCJL) 

Н2-Н1жж2 

H3=H1*H2 

H4=H2**2 

HS=H2*H3 

H6=H3**2 

мез. 

J2=J1+NM1 

J3=J2+NM1 

J4=J3+NM1 

YDX2-A(J1)/4.0*H4-A(J2) /3.0*H3-A(J3)/2.0*H2-A(J4) *H1 
BUG=-A(J1)/5.0*H5-A(J2)/4.0*H4-A(J3)/3.0*H3-A(J4)/2.0*H2 
XYDX=BUG+XIN(J1)*YDX 
BUG=-A(J1)/6.0*H6-A(J2)/5.0*H5-A(J3)/4.0*H4-A(J4)/3.0x*H3 
XXYDX=BUG+2.0*XIN(J1)*XYDX-XIN(J1)**2*YDX 

DO 7 N=JL,JU 

Hi=XIN(N+1)-XIN(N) 

IF(N.EQ.JU) Hi=XU-XIN(N) 

Н2=Н1**2 

H3=H1*H2 

H4=H2**2 

H5-H2*H3 

H6-H3**2 

Ji=N 

J2=J1+NM1 

J3=J2+NM1 

J4=J3+NM1 
BUG=A(J1)/4.0*H4+A(J2)/3.0*H3+A(J3)/2.0#H2+A(J4)*H1 
YDX=YDX+BUG 
CAT2A(J1)/5.0*H5*A(J2) /4.0*H4*A(J3) /3. 0*H3*A(J4) /2.0*H2 
PIG=CAT+XIN(J1)*BUG 

XYDX=XYDX+PIG 
DOG=A(J1)/6.0*H6+A(J2)/5.0*H5+A(J3)/4.0*H4+A(J4)/3.0+H3 
XXYDX=XXYDX+DOG+2.0*XIN(J1)*PIG-XIN(J1)**2*BUG 

CONTINUE 

RETURN 

END 


E а а-ы а-ы а-ы ран н 


С 


SUBROUTINE LINDK(NIN,XIN,YIN,A) 
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C----- GENERATES CUBIC COEFFICIENTS FOR PIECEWISE LINEAR FIT------------- 
REAL XIN(*),YIN(*) ,A(*) 
NM1-NIN-1 
DO i N=1,NMi 
A(N)=0.0 
M=N+NMi 
A(M)=0.0 
M=M+NM 1 
ACM)ZCYIN(N*1)-YIN (N)) /CXIN(N*1) -XIN(N) ) 
M=M+NM1 
1 A(M)=YIN(N) 
RETURN 
END 
ESCOCURGSGGOGCOSECCCCOSRGCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
SUBROUTINE PEVALDK(NIN,NOUT,ARK,S,XO, YO,CUB1,CUB2) 
REAL X0(200),Y0(200),S(200) , ARK(200) 
REAL CUBi((200-1)*4) ,CUB2( (200-1) #4) 
CALL EVALDK(NIN,NOUT,S,ARK,XO, CUB1) 
CALL EVALDK(NIN,NOUT,S,ARK, YO, CUB2) 
RETURN 
END 
We eGecececcc weccCceeeeecccccccccccccccccccccccccccccccccccccecccc 
SUBROUTINE PUGLYDK(NIN,X,Y,S,CUB1,CUB2) 
REAL X(200),Y(200),S(200) 
REAL CUB1((200-1)*4),CUB2((200-1) *4) 
STOFE0.0 
$(1)=0.0 
DO 10 I=2,NIN 
S(I)-SQRT((X(ID)-X(I-1)) *«24 (Y(I)-Y(I-1))**2)-STOT 
STOT=STOT+SQRT((X(CI)-X(I-1))**2+(Y(I)-Y(I-1))**2) 
10 CONTINUE 
DO 20 I=2 NIN 
S(I)=S(I)/STOT 
20 CONTINUE 
NCL=1 
NCR=1 
ESL=0 
ESR=0 
CALL UGLYDK(NIN,NCL,NCR,S,X,ESL,ESR,CUB1) 
CALL UGLYDK(NIN,NCL,NCR,S,Y,ESL,ESR,CUB2) 
RETURN 
END 


SUBROUTINE UGLYDK(NIN,NCL,NCR,XIN,YIN,ESL, ESR, AE) 
en 1975 DUCK SERIES J.E.KERWIN MODIFIED 6/21/82--------------------- 


REAL XIN(*),YIN(*) ,AE(*) ,H(200),D(200) ,AU(200) , AM(200), 
^ S(200) ,AL(200),X(200) 

DATA HALF/0.5E00/,TW0/2.0E00/,SIX/6.0E00/,RAD/1.745329E-02/ 
NM1=NIN-1 

NM2=NM1-1 

NM3=NM2-1 

NEQ=NM2 

DO 1 N=1,NM1 

H(N)=XIN(N+1)-XIN(N) 

1 D(N)=(YIN(N+1)-YIN(N))/H(N) 

IF(NCL.EQ.2) NEQ=NEQ+1 

IF(NCR.EQ.2) NEQ=NEQ+1 

NSQ=NEQ+**2 

4-1 

IF(NCL.LT.2) GO TO 6 

AM(1)=TWO*H(1) 

AU(1)=H(1) 

SLP-ESL*RAD 

S(1)=(D(1)-TAN(SLP))*SIX 
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22 


10 


J=J+1 
AL(2)=H(1) 

DO 5 N=1,NM2 

IF(N.GT.1) AU(J-1)=H(N) 

AM(J)=TWO*(H(N)+H(N+1)) 

IF(N.LT.NM2) AL(J+1)=H(N+1) 
IF(N.EQ.2. AND. NCL. EQ. 1) AU(J-1)-AU(J-1) -H(N-1) **2/H(N) 
IF(N.EQ.1.AND.NCL.EQ.1) AMC(J)=AM(J)+(1.0+H(N)/H(N+1))*H(N) 
IF(N.EQ.NM2.AND.NCR.EQ.1) AM(J)=AM(J)+(1.0+H(N+1)/H(N))*H(N+1) 
IF(N.EQ.NM3.AND.NCR.EQ.1) AL(J+1)=AL(J+1)-H(N+2)**2/H(N+1) 
S(I)=(D(N+1)-D(N))*SIX 

J=J+1 

CONTINUE 

IF(NCR.LT.2) GO TO 7 

AL(NEQ)=-H(NM1) 

AM(NEQ)=-TWO*H(NM1) 

AU(NEQ-1)=H(NM1) 

SLP-ESR*RAD 

S(J)=(D(NM1)+TAN(SLP))*SIX 

CONTINUE 

DO 4 K=2,NEQ 

AL(K)=AL(K)/AM(K-1) 

AM(K)=AM(K)-AL(K)*AU(K-1) 
S(K)=S(K)-AL(K)*S(K-1) 

CONTINUE 

X (NEQ)=S(NEQ)/AM(NEQ) 

DO 2 L=2,NEQ 

K=NEQ-L+1 

X(K)=(S(K)-AU(K)*X(K+1))/AM(K) 

CONTINUE 

DO 22 N-1,NEQ 

S(N)=X(N) 

HOLD=S(NEQ) 

IF(NCL.EQ.2) GO TO 8 

DO 9 N=1,NM2 

M=NM2-N+2 

S(M)=S(M-1) 

IF(NCL.EQ.0) S(1)=0.0 

BUG=H(1)/H(2) 

IF(NCL.EQ.1) S(1)=(1.0+BUG)*S(2)-BUG*S(3) 
IF(NCR.EQ.0) S(NIN)=0.0 

BUG=H(NM1)/H(NM2) 

IF(NCR.EQ.1) S(NIN)=(1.0+BUG)*S(NM1)-BUG*S(NM2) 
IF(NCR.EQ.2) S(NIN)=HOLD 

DO 10 N=1,NM1 

AE(N)=(S(N+1)-S(N))/(SIX*H(N)) 

M=N+NM1 

AE(M)=HALF*S(N) 

M=M+NM1 

AE(M)=D(N)-H(N) *(TWO*S(N)+S(N+1)) /SIX 

M=M+NM1 

AE(M)=YIN(N) 

RETURN 

END 
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Appendix B 


sample FIT2D Input Files 


B.1 Bounded Foil (fname.ctrl) 


Header: Sample for Cupped B-1 

Foil geometry file name 

reup foil 

AOA Xpiv Ypiv 

poc 0.3 0.01 

Chord/Tunnel Width 

Duo 

USE DSL PHI1 РНІ2 

2-279 .0.5 0.25 

NUS NDS NVS NTOP NBOT(Always 8*n-1 points) 

Dn. 0055 95 95 

RESLE RESMID RESTE RESWAL  RESBL PACK 

0.001 0.03 0.001 1.0e-5 1.0е-5 0.20 

Re#(chord based) 

3e6 

Desired Inmesh Input file: 

cupi.dat 

Desired Inmesh Restart file 

cupr.dat 

Number of INMESH Iterations: 

1500 

Convergence tolerance: 

1.0e-10 

Wake Data(-1 none, O VLM, or nranswk) 
10 

.713324 -0.016042 

28317 -0.022880 

.819395 -0.032600 

.027451 -0.043110 

. 210484 -0.052885 

.424693 -0.061495 

.674272 -0.069320 

„353325 -0.077086 

72273109 -0 083998 

-580143 -0.089997 


мМ м нъ њан Оо Оо о 
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B.2  Unbounded Foil (fname.ctrl) 


Header: Sample Unbounded Cupped B-1 Foil 

Foil geometry filename: 

reup.foil 

ADA Xpıv Ypiv 

што. 3 0.01 

Chord/Tunnel Width 

0.08 

ШЕШЕ ШЕТ, РНІ1 РНҮ2 

EEOE5.0. 2.00 1.00 

NUS NDS NVS NTOP NBOT(Always 8*n-1 points) 

er er 95 95 95 

RESLE RESMID RESTE RESWAL  RESBL PACK 

0-001 0.03 0.001 0.005 1.0e-5 0.20 

Re#(chord based) 

3e6 

Desired Inmesh Input file: 

cupi.dat 

Desired Inmesh Restart file 

cupr.dat 

Number of INMESH Iterations: 

1500 

Convergence tolerance: 

1.0e-10 

Wake Data(-1 none, O VLM, or nranswk) 
20 

.713324 -0.016042 

- АТО .022880 

309335720 .032600 

.027451 -0.043110 

.210484 -0.052885 

.424693 -0.061495 

.674272 -0.069320 

2953225 -0.077086 

2251516 -0.083998 

.580143 -0.089997 

ЭКМЕШЕВЕ-0. (095028 

. 201140 -0 099227 

.603834 -0.102689 

. 939133 -0.105482 

261238 -0.107661 

1565296 -0.109321 

.846685 -0.110522 

.101448 -0.111331 

.326634 -0.111836 

lease 0.112118 


сл сла н> ОО О) 0 М мМ н ҥн к н о Оо О 
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B.3 Sample Foil Geometry File (fname.foil) 


Offsets start at the trailing edge marching forward on lowersurface, around the leading 
edge and along upper surface to trailing edge. It is not neccessary to have an explicit 
point at the leading edge. FIT2D splines the offsets and finds the leading edge. 


FIT2D assumes that this input file is at zero degrees angle of attack. 


Cupped B-1 Foil 
2 


720 

161 

1.000000 0.000000 
0.988517 0.002052 
022791109 0.003558 
0.968851 0.004586 
0.958810 0.005204 
0.949059 0.005478 
0.939668 0.005477 
0.930708 0.005269 
0.922250 0.004920 
0.914365 0.004499 
0.007665 -0.006030 
090-590 -0.005229 
0.003652 -0.004320 
0.002127 -0.003289 
0.000963 -0.002132 
0.000220 -0.000846 
0.000041 0.000570 
0.000133 0.001781 
0.000676 0.003074 
0.001617 0.004450 
0.002986 0.005909 
0.004814 0.007453 
0.007129 0.009082 
0.009962 0.010796 
029352072 0.028109 
0.962481 0.024014 
0.971452 0.019244 
0.980651 0.013711 
0.990145 0.007326 
1.000000 0.000000 
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Appendix C 


Marine Hydrodynamics Lab Water 
Tunnel Geometry 


C.1 System Overview 





Figure C-1: MIT MHL Water Tunnel 
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Figure C-2: MIT MHL Water Tunnel Test Section 
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Appendix D 


GETWAKE Program Listing 


Q000€000050500000 


Q 


оо 


98 


PROGRAM GETWAKE 


Reads in a set of wakepoints . 


INPUT REQUIREMENTS: 

JUNK header line 

NPOINTS - NUMBER OF OFFSETS 
XIN, YIN 


(to end of file) 


IMPLICIT NONE 
REAL XIN(1000),YIN(1000) 
REAL Xout(200),Yout (200) 
integer npoints,nfoil,i,k,ntemp,nout,j,ierr 
character FNOPEN*20 
character junk*30, title*30, JUNK2*2 
character PROMPT2*30 
character FIN*20 
WRITE(*,*) ?*xx? 
write(*,’(A,$)’) ’Enter Filename of Tecplot stream data: ’ 
read(*,’(A)’) FNOPEN 
WRITE(*,*)?xx*x? 
write(*,’(A)’) ’Shortening file:’ // FNOPEN 
OPEN (UNIT=1 ,FILE=FNOPEN , FORM=’ FORMATTED’ , STATUS=’OLD’ ) 
WRITE(*,*)? xxx? 
WRITE(*,*)? ***READING DATA ІМжжж? 
WRITE(*,*)?xxx? 
read(1,'(A)?) junk 
read(1,'(A)?) junk 
read(1,’(A)’) junk 
WRITE(*,*)'READ THE TOP OF THE FILE? 


keep next read as the title 
read(1,’(A)’) title 


now want to extract the integer for number 
of pairs of points 
CALL READIJ(1,npoints,J,IERR) 
FORMAT (A23,16,A5S) 
WRITE(*,98)’***NUMBER DATAPOINTS IN: ’,NPOINTS,’ ***? 


WRITE(*,*)?xxx? 
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100 


91 
92 


200 


READ(1,’(A)’) JUNK 
READ(1,*) ((xin(k),yin(k)),K=1,npoints) 
READ(1,*) JUNK 
CLOSE(1) 
WRITE(*,*(A)?)” ***CLOSING FILE: * // FNOPEN // ? жжж› 
WRITE(*,*x)? xxx? 


Shorten Data file by taking only 20 datapoints in range 


NOUT=20 
WRITE(*,98)'***NUMBER OF POINTS OUT: ?,NOUT,? xxx? 
WRITE(*,*x)?xxx? 


NTEMP-int(float(npoints)/NOUT) 
WRITE(*,98)?'***DATA SKIP INTERVAL: ',NTEMP,? xxx? 
WRITE(*,*) ’ xxx? 


do 100 i=1,NOUT 
xout(i)=xin(i*ntemp-ntemp+1) 
yout(i)=yin(i*ntemp-ntemp+1) 
continue 


write(*,*)’****Data has been shortened.*xx*xxx? 
Write(*,*)?xxx? 


PROMPT2 = ’Output FILE’//’ (’//’wakelin.dat’//’) =? 
WRITE (*,’(A,$)’) PROMPT2 

READ (*,’(A)’) FIN 

IF (FIN(1:1).EQ.’ ’) FIN = ’wakelin.dat’ 

WRITE (*,’(A)’) ’OPENING FILE: ^ // FIN 

OPEN (UNIT23,FILE-FIN,STATUS- UNKNOWN ›) 


format(i4) 
format(2f10.6) 
WRITE(x*,*)?xxx? 
WRITE(*,*)'***WRITING OUTPUT TO FILEx*x' 
WRITE(*,*x)?xxx? 
WRITE(3,*)’Shortened: ’,TITLE 
WRITE(3,91)NOUT 
do 200 i=1,nout 

write(3,92) xout(i),yout(i) 
continue 
close(3) 
write(*,*)’ All Done, Thanks. ? 
WRITE(*,*) xxx? 

stop 

end 


CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 


с 
с 


жжжжжжжжжжжжжжжжжжжжжжжжжжж SUBROUTINE READIJ жжжжжжжжжжжжжжжжжжжжж 
EXTRACTS THE INTEGER I,J INDICES FROM A TECPLOT FILE HEADER * 


ж 


ж 
к 
к 
к 
ж 
ж 
ж 


Arguments: IUNIT....Unit number of Tecplot file 


о Number of columns in IJ ordered data 
er Number of rows in IJ ordered data 
TERR..... O=Valid data, 1=Subroutine failed 


It is assumed that I and J are between 1 and 9999 
Justin E. Kerwin August 19,1993 


ж 


Ж 
ж 
Ж 
* 
* 


ЖЖЖЖЖЖЖЖ KK KH EEK FF FF FF FE EEE € 


SUBROUTINE READIJ(IUNIT,I,J,IERR) 
CHARACTER*80 LABEL 

CHARACTER*4 ICODE 

ICODE=’ ) 
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IERR=O 

E. Read in the Tecplot header line as a character string LABEL. 
READ(IUNIT, '(A)?) LABEL 

es Endaenedlenseth of the character string.........-. s. rrr 
LMAX-LEN(LABEL) 
LBEGIN-1 

Bo БИЕСІН 1: (1J=1), then ҒІШЕ 7:(17-2)..................... 


DO 


100 


э ә ә ә е ее o 


200 


300 


400 IJ-1,2 
ICODE=’ | 
Find the position of the first = sign in the string..... 
DO 100 L=LBEGIN,LMAX 
IF(LABEL(L:L).EQ.CHAR(61)) THEN 
LMIN=L+1 
GO TO 110 
END IF 
CONTINUE 
IERR=1 
GO TO 9999 
CONTINUE 
Find the position of the first number following an = sign. 
DO 200 L=LMIN,LMAX 
IF(LABEL(L:L).GT.CHAR(48).AND.LABEL(L:L).LT.CHAR(58)) THEN 
LSTART=L 
GO TO 210 


GO TO 9999 
K=1 
Generate a substring ICODE consisting of consecutive numbers 
DO 300 L=LSTART,LSTART+4 
IF(LABEL(L:L).LT.CHAR(48).OR.LABEL(L:L).GT.CHAR(57)) THEN 
LBEGIN=L 
GO TO 310 
ELSE 
ICODE(K:K)=LABEL(L:L) 
K=K+1 
END IF 
CONTINUE 
IERR=1 
GO TO 9999 
Convert the substring to integers I,J using an internal read 
IF(IJ.EQ.1) READ(ICODE,'(I4)?) I 
IF(IJ.EQ.2) READ(ICODE,’(I4)’) J 


400 CONTINUE 


9999 RETURN 


END 
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Appendix E 


PATCH Progam Listing 


PROGRAM PATCH 
REAL*8 X(200,200),Y(200,200) 
character*10 FNOPEN 
FNOPEN-'inpO2.dat' 
OPEN (UNIT=1 , FILE=FNOPEN, FORM=’ FORMATTED’ , STATUS=’ OLD’ ) 
READ(1,*) JGRD,KGRD 
READ(1,*) ((X(J,K),J=1,JGRD) ,K=1,KGRD) 
READ(1,*) ((Y(J,K),J=1,JGRD),K=1,KGRD) 
CLOSE (1) 


FNOPEN=’inpO2.dat’ 
OPEN(UNIT=2,FILE=FNOPEN, FORM=’FORMATTED’ ‚STATUS=’UNKNOWN’ ) 
WRITE(2,*) JGRD-1,KGRD 

WRITE(2,*) ((X(J,K),J=2,JGRD),K=1,KGRD) 

WRITE(2,*) ((Y(J,K),J=2,JGRD),K=1,KGRD) 

CLOSE(2) 


FNOPEN=’ inp03.dat’ 

OPEN (UNIT=1 , FILE=FNOPEN , FORM=’ FORMATTED’ , STATUS=’ OLD’ ) 
READ(1,*) JGRD,KGRD 

READ(1,*) ((X(J,K),J=1,JGRD),K=1,KGRD) 

READ(1,*) ((Y(J,K),J=1,JGRD),K=1,KGRD) 

CLOSE (1) 


FNOPEN-' inp03.dat' 
OPEN(UNIT-2,FILE-FNOPEN,FORM-FORMATTED' , STATUS-' UNKNOWN ' ) 
WRITE(2,*) JGRD-1,KGRD 

WRITE(2,*) ((X(J,K),J=2,JGRD),K=1,KGRD) 

WRITE(2,*) ((Y(J,K),J=2,JGRD),K=1,KGRD) 

CLOSE(2) 


FNOPEN=*inp05.dat? 

OPEN (UNIT=1, FILE=FNOPEN , FORM=’ FORMATTED’ , STATUS=’ OLD’ ) 
READ(1,*) JGRD,KGRD 

READ(1,*) ((X(J,K),J=1,JGRD),K=1,KGRD) 

READ(1,*) ((Y(J,K),J=1,JGRD),K=1,KGRD) 

CLOSE(1) 


FNOPEN-'inpO5.dat' 
OPEN(UNIT-2,FILE-FNOPEN,FORM-?FORMATTED ' , STATUS- ' UNKNOWN? ) 
WRITE(2,*) JGRD-1,KGRD 

WRITE(2,*) ((X(J,K) ,J=2,JGRD) ,K=1,KGRD) 

WRITE(2,*) ((Y(J,K),J=2,JGRD),K=1,KGRD) 

CLOSE (2) 
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FNOPEN=’ inp06.dat’ 
OPEN(UNIT=1, FILE=FNOPEN, FORM=’FORMATTED’ , STATUS=’OLD’ ) 
READ(1,*) JGRD,KGRD 

READ(1,*) ((X(J,K),J=1,JGRD),K=1,KGRD) 

READ(1,*) ((Y(J,K),J=1,JGRD),K=1,KGRD) 

CLOSE(1) 


FNOPEN-' inp06.dat' 
OPEN(UNIT=2,FILE=FNOPEN,FORM=*FORMATTED? ,STATUS='UNKNOWN? ) 
WRITE(2,*) JGRD-1,KGRD 

WRITE(2,*) ((X(J,K),J22,JGRD) ,K21,KGRD) 

WRITE(2,*) ((Y(J,K),J22,JGRD) , KZ1,KGRD) 

CLOSE(2) 


stop 
end 
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Appendix F 


Case Study Foil Offsets 


Offsets start at the trailing edge of the pressure side and march forward around the 
leading edge to the trailing edge of the suction side. Both data sets use a normalized 
chord length 1.0. The offsets listed here correspond to the foils displayed in Figures 


5-1 and 5-9. 
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КЛ 


ооо ооо ооо оо оо оо ооо оо оо оо ооо оо оо о ооо, оо о ооо, оо, о, ооо оо, 


Foil 


20 0 
. 9245» —0 
23907209 
.985 -0 
.975 -0 
1965 20. 
о c 
2925 -0 
.900 0. 
.875 0 
2850 -0 
4525 -0% 
.800 -0 
LORO. 
790 =0. 
„125 -0. 
100 =0. 
650 =0. 
.600 -0. 
.550 -0. 
.500 -0. 
.450 -0. 
.400 -0. 
2529 -0 
.300 -0. 
.290 -0. 
„202 -0. 
„15 -0. 
:150 -0. 
„123 -0. 
-100 -0. 
ОТЫ =0. 
.050 -0. 
1035 =0. 
-025 =0. 
1015 =0 
-010 =0 
.005 -0 
.0025 -0.003274 
2001420 
.0 0 
.001 0.002076 
.0025 0 
.005 
.010 
‚015 


I 
ол 
о 
оэоэооооооооооооосоо« 


vO 


.000896 
.000797 
.000703 
.000537 
000397 
.000225 
.000031 
000049 


.000009 


.000150 
000405 
.000734 
001126 
001574 
002076 
002631 
003968 
005685 
007780 
010205 
012882 
015697 
.018500 
021093 
023204 
024436 
024512 
024077 
023054 
021373 
018969 
015681 
013111 
011017 
.008420 
.006793 
.004709 


.002036 


.0 


.003375 


.004911 
.007200 
.009036 
2042061 
.014599 
.017865 
. 022406 
.026204 
.029474 
10352329 
.034841 
.037059 
.040736 
.043546 
.045599 
.046962 
.047677 
.047760 


Case Study I: The HRA 
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г“оооооооооосоосооооооо 


ооо ооо оо ооо осо 


. 047207 
. 045990 
.044037 
. 041199 
.039327 
.037072 
. 034427 
.031406 
. 028036 
. 024371 
.020515 
‚016571 
‚012588 
. 008599 
.006225 
.004671 
.003166 
.002433 
1001713 
„О 





F.2  B-1 Foil With Cup Mod- 
ification To Trailing Edge 


OOOOOOOO0O0O0O0O0000000000000000000000000o0ooooooooooooooooooooooon 


. 000000 
. 989547 
.979110 
.968851 
.958810 
.949059 
.939668 
.930708 
. 922250 
.914365 
3907 1294 
.900497 
.894374 
.888628 
2883131 
‚877761 
‚872391 
. 866896 
2991152 
.855033 
.848423 
.841250 
.833456 
.824986 
. 815788 
„805789 
‚794947 
193201 
.770495 
.756770 
.742000 
. 726196 
. 709371 
. 691542 
2672721 
. 652924 
. 632166 
.610461 
.587823 
.564277 
.539906 
.514822 
.489138 
.462964 
.436413 
.409596 
7962625 
-355611 
‚328007 
2304935 
‚275618 
‚249922 
. 225056 
201226 
.178641 
.157509 
. 138036 
. 120430 
. 104875 
1091313 
. 079551 


ооо ооо оо о ооо оооо 


. 000000 
.002052 
.003558 
.004586 
. 005204 
.005478 
. 005477 
.005269 
.004920 
. 004499 
. 004072 
.003672 
. 003298 
.002943 
.002604 
. 002274 
. 001950 
.001626 
. 001297 
. 000959 
. 000606 
. 000236 
. 000154 
.000568 
. 001006 
.001473 
.001971 
.002503 
.003071 
.003678 
.004322 
.004997 
.005694 
. 006406 
.007126 
.007845 
.008556 
2009251 
.009922 
.010563 
.011168 
.011734 
‚012255 
‚012730 
013153 
2013521 
.013830 
.014076 
.014255 
. 014366 
. 014414 
.014402 
. 014337 
.014221 
.014061 
2913801 
‚013626 
„013359 
‚013067 
012752 
.012417 
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OOOOOO0O0O0O0O0000000000000000000000000000000000000o0o0ooooooooooooooooo 


.069396 
.060654 
2053133 
.046637 
.040974 
.035949 
.031370 
.027095 
.023110 
.019418 
.016024 
012931 
.010143 
.007665 
. 005500 
. 003652 
. 002127 
. 000963 
.000220 
. 000000 
. 000133 
. 000676 
. 001617 
.002986 
. 004814 
. 007129 
. 009962 
. 013344 
2017305 
.021881 
.027105 
2933011 
.039635 
.047010 
.055170 
.064149 
- 073983 
. 084704 
. 096348 
. 108948 
. 122492 
130929 
.152161 
.168152 
. 184827 
2202 DI 
.219958 
249209 
.257024 
„276117 
‚ 295494 
-315089 
. 334845 
. 354719 
. 374667 
. 394646 
.414616 
.434532 
.454353 
.474035 
.493537 
‚512816 
‚531825 
2550533 
.568899 


.012064 
.011696 
011315 
.010923 
.010524 
.010120 
.0097 13 
.009301 
.008873 
.008416 
.007916 
. 007361 
.006736 
. 006030 
.005229 
.004320 
.003289 
. 002132 
.000846 
. 000570 
.001781 
.003074 
.004450 
.005909 
.007453 
.009082 
.010796 
.012596 
.014480 
.016443 
.018482 
.020591 
.022766 
.025003 
2022297 
. 029643 
.032038 
.034476 
. 036953 
.039464 
. 042000 
.044547 
.047087 
. 049605 
. 052085 
.054513 
.056871 
.059145 
‚061319 
‚063376 
‚065302 
. 067080 
. 068702 
.070169 
.071483 
.072647 
.073662 
.074531 
.075256 
. 075838 
.076282 
.076587 
.076757 
.076795 
.076701 





"ЭӘОЭЭООООО ООЭЭООО ОО ОО ООООО ОВ ОоОооооооооо 


. 986921 
.604596 
.621920 
.638890 
.655504 
.671758 
.687648 
-7/03172 
- ( 18827 
. 733108 
.747514 
.761542 
„15192 
. 188463 
.801358 
‚813876 
‚826017 
79527162 
.849171 
.860186 
.870825 
.881089 
.890979 
.900510 
.909739 
„918732 
921555 
. 936276 
1944959 
ISSO TZ 
.962481 
. 971452 
.980651 
990145 
.000000 


оээЭэОЭОооооооооооооооооооооооооооооооо 


ӨЛСЕ 7 
‚076125 
‚075647 
‚075044 
‚074317 
‚073469 
‚072500 
‚071413 
. 070208 
.068888 
.067454 
.065910 
.064268 
‚062543 
‚060746 
‚058891 
‚056991 
‚055059 
‚053108 
„051152 
‚049203 
‚047275 
‚045380 
‚043511 
. 041590 
.039530 
‚037241 
. 034634 
.031619 
2028109 
-024014 
20193272 
.013711 
.007326 
. 000000 
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